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ABSTRAC 


This  report  gives  a  summary  of  the  projects  conducted 
during  the  last  year  under  the  contract  Research  in  Seis¬ 
mology.  These  investigations  fall  roughly  under  three 
broad  headings:  (1)  Theoretical  studies  of  source  mech¬ 
anisms  of  earthquakes  and  underground  nuclear  explosions; 

(2)  Earth  structure  and  path  effects,  especially  the  effects 
of  lateral  heterogeneities  in  velocity  and  in  Q;  (3)  Studies 
of  the  capabilities  of  arrays  for  event  detection  and  loca¬ 
tion.  Lists  of  publications  and  theses  completed  during  the 
•contract  year  are  also  included. 
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SUMMARY 


In  this  annual  report  we  review  the  work  completed 
under  the  contract  Research  in  Seismology  during  the  year 
1  April  1972  through  31  March  1973.  Within  the  broad  guide¬ 
lines  of  the  problem  of  discriminating  earthquakes  from 
underground  nuclear  explosions,  we  have  conducted  a  number 
of  specific  investigations  of  seismic  sources  and  seismic- 
wave  propagation. 

The  topics  studied  can  be  grouped  under  three  broad 
headings : 

(1)  Source  mechanisms  of  earthquakes  and  explosions, 

(2)  Earth  structure  and  path  effects, 

(3)  Array  studies. 

In  the  following  sections  we  present  abstracts  of  papers 
published  or  soon  to  be  published  that  fall  into  each 
category.  Recently  completed  work  is  discussed  in  fuller 
detail.  We  also  include  lists  of  all  publications  and 
theses  supported  under  this  project  during  the  contract  year. 

Several  theoretical  approaches  to  seismic  source  mech¬ 
anisms  have  been  followed.  Y.  Ida  has  pursued  further  his 
propagating-crack  model  for  earthquake  rupture.  The  model 
has  the  advantage  that  it  links  seismic  observations  such 
as  near-field  acceleration  to  physical  properties  of  rocks 
such  as  static  and  dynamic  friction.  Variations  in  the 
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parameters  used  to  formulate  the  boundary  condition  across 
the  fault  plane  can  explain  the  transition  between  stick- 
slip  and  stable  sliding  observed  in  rock  failure,  the  dis¬ 
similarity  in  the  amplitude  spectra  of  large  and  small 
earthquakes,  and  rates  of  seismicity  along  faults. 

A  number  of  related  theoretical  seismic  source  models 
were  constructed  by  S.  Singh  His  approach  was  to  super¬ 
pose  point  forces  along  some  boundary  in  a  medium  of  speci¬ 
fied  elastic  or  viscoelastic  properties  and  to  solve  for 
the  resulting  displacement  and  stress  fields.  Among  the 
interesting  results  of  his  studies  are  the  demonstration 
that  apparently  torsion-free  sources  can  generate  SH-type 
waves  and  that  the  effects  of  viscoelasticity  can  signifi¬ 
cantly  alter  the  displacement  of  stress  fields  for  certain 
kinds  of  faulting  from  those  calculated  assuming  perfect 
elasticity. 

A  number  of  advances  were  made  toward  inverting  normal¬ 
mode  and  travel-time  data  to  infer  models  of  earth  structure. 
Bounds  on  shear  velocity  and  density  in  the  mantle  and  core 
were  obtained  explicitly  by  C.  Johnson  using  a  new  linear 
programming  technique.  In  a  related  work,  Wiggins,  McMechan 
and  Toksflz  developed  a  procedure,  based  on  the  Herglotz- 
Wiechert  integral,  for  direct  determination  of  the  envelope 
of  possible  body-wave  velocities  in  the  mantle  and  core. 
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Both  of  these  studies  have  both  clarified  the  limits  of 
uncertainty  in  current  earth  models  and  pointed  toward  new 
data  sets  that  would  improve  present  constraints. 

A  major  focus  of  our  attention  has  been  directed  toward 
understanding  lateral  as  well  as  radial  variations  in  earth 
structure.  As  a  consequence  of  a  P-wave  travel-time  study 
using  waves  from  deep  earthquakes  to  minimize  near-source 
heterogeneity,  M.  Sengupta  and  B.  Julian  noticed  systematic 
differences  in  deep-mantle  travel  times  for  various  paths. 

A  lateral  variation  of  P-wave  velocity  in  the  lower  mantle 
of  at  least  1  percent  is  required. 

The  anomalous  patterns  of  travel-time  delays  and  ampli¬ 
tude  variations  across  LASA  have  been  recognized  for  some 
time.  K.  Aki  has  shown  that  this  may  be  attributed  to 
scattering  by  random  velocity  fluctuations  in  the  crust 
beneath  the  array.  The  scattering  is  quite  severe  and  the 
simple  Chernov  theory  which  fits  the  observations  at  lower 
frequencies  (.5  Hz)  is  inadequate  at  higher  frequencies 
(>  1  Hz)  . 

The  theory  of  plate  tectonics  has  led  to  the  recognition 
of  a  number  of  specific  heterogeneities  in  earth  structure, 
in  particular  the  regions  in  the  vicinity  of  plate  boundaries 
at  spreading  centers  and  at  subduction  zones .  When  seismic 
sources  are  located  in  such  regions,  the  effects  on  the 
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propagation  of  body  waves  can  be  pronounced.  S.  Solomon 
has  demonstrated  the  striking  azimuthal  dependence  of  shear- 
wave  spectra  from  earthquakes  near  ridges  due  to  a  localized 
zone  of  very  low  Q  in  the  uppermost  mantle  beneath  the 
ridge  crest.  M.N.  Toksflz  and  his  colleagues  N.  Sleep  and 
A.  Smith  have  further  pursued  their  models  of  temperature 
and  stress  in  subducted  slabs  of  lithosphere.  Current 
models  successfully  explain  the  travel-time  and  amplitude 
anomalies  associated  with  shallow  earthquakes  and  explosions 
on  island  arcs  and  the  distribution  and  focal  mechanisms 
of  deep,  intra-slab  earthquakes. 

The  efficient  utilization  of  the  large  aperture  seis¬ 
mic  arrays  LASA  and  NORSAR  has  long  been  one  of  our  continuing 
interests.  In  his  Ph.D.  thesis  and  later  work,  S.  Shlien 
has  made  extensive  statistical  studies  of  the  capability  of 
the  combined  arrays  for  on-line  signal  detection  and  event 


location . 
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2.  SOURCE  MECHANISMS  OF  EARTHQUAKES  AND  EXPLOSIONS 

2.1.  Stress  Concentration  and  Unsteady  Propagation  of  Longi¬ 
tudinal  Shear  Cracks  by  Y.  Ida  (Abstract) 

Unsteady  propagation  of  a  longitudinal-shear  crack  is 
studied  through  numerical  calculation  for  better  understand¬ 
ing  of  fracture  and  earthquakes.  The  relation  which  formu¬ 
lates  rupture  velocity  as  a  function  of  stress  intensity  is 
used  to  describe  the  fracture  process  at  the  crack  tip.  The 
effects  of  creep  and  frictional  slidings  are  taken  into 
account  by  assuming  a  suitable  form  of  the  slip-velocity 
dependent  boundary  condition.  In  this  model,  stress  accumu¬ 
lation  takes  place  through  creep,  and  then  the  crack  begins 
to  grow,  controlled  by  the  dynamic  friction  across  the  crack 
surface.  It  is  shown  that  the  motion  of  the  crack  tip  is 
either  smooth  or  bumpy,  depending  on  the  parameters  of  the 
boundary  condition.  The  brittleness  and  ductility  of  mate¬ 
rial  are  interpreted  by  these  two  types  of  crack  propagation, 
and  this  interpretation  is  qualitatively  consistent  with  the 
experimental  results  on  brittle-ductile  transition  induced 
by  the  temperature,  the  confining  pressure,  and  the  pore 
pressure.  It  is  shown  that  the  presence  of  smooth  and  bumpy 
crack  propagations  might  give  a  possible  explanation  on  the 
dissimilarity  between  large  and  small  earthquakes,  pointed 


out  by  Aki.  This  model  also  suggests  a  method  of  estimating 
the  recurrence  period  of  earthquakes,  based  on  the  crustal 
deformation  across  the  fault. 

2 . 2  Cohesive  Force  and  Unsteady  Propagation  of  a  Longi¬ 
tudinal-Shear  Crack  by  Y .  Ida  (abstract) 

The  dynamic  propagation  of  a  longitudinal-shear  crack 
is  studied  to  form  a  model  of  fracture  that  has  some  physi¬ 
cal  basis.  The  stress  field  around  the  crack  tip  is  analyzed 
in  detail,  .der  the  assumption  that  the  cohesive  force  is 
given  as  a  function  of  the  displacement  discontinuity.  This 
analysis  yields  a  fracture  condition  which  turns  out  to  be 
equivalent  to  the  energetic  criterion,  but  the  meaning  of 
specific  surface  energy  is  made  clearer.  The  stress  distri¬ 
bution  in  the  vicinity  of  the  crack  tip  is  demonstrated  for 
several  models  of  cohesive  force  diagrams.  It  is  assumed 
that  the  slip  velocity-dependent  boundary  condition,  which 
includes  the  effect  of  creep  and  the  static  and  dynamic 
frictions  governs  the  semi-infinite  fault  plane  except  for 
the  infinitesimally  small  end-region  associated  with  cohe¬ 
sive  force.  Based  on  these  formulations  of  the  fracture 
and  boundary  conditions,  the  crack  growth  is  calculated 
numerically.  In  this  model,  the  stress  accumulation  first 
ikes  place  through  the  creep  over  the  whole  fault  plane. 
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and  then  the  crack  begins  to  grow,  accompanied  by  the  for¬ 
mation  of  the  crack  surface  havinq  the  dynamic  friction. 

It  is  shown  that  the  crack  tip  is  either  "smooth"  or 
"bumpy,"  depending  on  the  parameter  of  the  boundary  con¬ 
dition.  Smooth  propagation  is  associated  with  the  rapid 
expansion  of  the  crack  surface,  while  the  bumpy  propagation 
only  produces  the  crack  surface  on  the  restricted  area 
adjacent  to  the  crack  tip.  The  brittle  and  ductile  frac¬ 
tures  of  the  material  are  interpreted  by  these  two  types 
of  crack  propagation,  and  this  interpretation  successfully 
predicts  the  brittle-ductile  transition  induced  by  the 
temperature,  the  confining  pressure,  and  the  pore  pressure. 

9 

2 . 3  The  Maximum  Acceleration  of  Seismic  Ground  Motion 
by  Y.  Ida 

Abstract 

The  near-field  particle  motion  is  evaluated  on  the 
basis  of  the  idea  that  the  initial  shape  of  the  seismic 
source  function  should  be  determined  by  the  fracturing 
process  in  the  rupture  front.  The  effect  of  the  cohesive 
force  on  the  source  function  is  examined  in  detail,  and  the 
following  relations  are  derived  for  the  maximum  velocity 
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nM  and  the  maximum  acceleration  u„;  U,.  ~  (o  /u)c  and 

M  M  M  O 

2  2 

~  (°0A‘)  (c  /R0^  '  where  oq  is  the  strength,  Pq  is 
the  displacement  discontinuity  required  for  the  fracture, 
c  is  the  rupture  velocity  and  p  is  the  rigidity.  If 
the  strength  of  rock  is  assumed  to  govern  the  earthquake 
rupture  (oq  '  1  kbar  and  Dq  -  10  cm)  ,  reasonable  value:, 
are  obtained,  as  UM  ~  100  ern/s  and  ~  1  g.  Th<.  peri  I 

and  the  dimension  associated  with  the  large  par  tic  Jo  acc»  J  •• 
eration  are  roughly  0.1  second  and  300  meters,  respectively, 
urder  the  above  assumption.  This  result  suggests  (he  possi¬ 
bility  that  eartaqu.ikc:  hazard  may  be  predicted  from  the  mech¬ 
anical  property  rocks. 
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INTROUUCTION 

In  the  vicinity  of  a  seismic  fault,  it  is  expected 
that  the  ground  motion  is  primarily  governed  by  the  fault- 
slip  time  function  rather  than  the  geometry  and  size  of 
the  fault  plane.  For  the  purpose  of  engineering  seismology, 
it  is  thus  very  important  to  determine  tl.o  accurate  form 
of  the  fault-slip  time  function.  Brune  (1970)  proposed  a 
simple  model  of  seismic  source,  assuming  that  the  slip 
motion  occurs  instantaneously  over  the  whole  fault  plane. 

In  his  model,  the  particle  motion  is  related  to  the  pre¬ 
existing  stress  field.  It  is  physically  more  natural, 
however,  to  regard  the  earthquake  occurrence  as  spontaneous 
rupture  propagation  associated  with  1  he  stress  concentration 
in  the  rupture  front.  If  this  is  the  case,  the  rupture 
cannot  propagate  supersonically,  and  Brune' s  time  function 
is  no  longer  applicable  (Ida  and  Aki,  1972). 

In  the  case  of  a  subsonic  rupture  propagation,  the 
entire  form  of  source  function  involves  several  complicated 
effects,  such  as  the  initial  state  of  the  displacement 
field,  and  the  dynamic  processes  of  fracture  and  friction. 
If  we  restrict  our  consideration  to  the  initial  rise  of 
source  time  function,  however,  the  analysis  can  be  made 
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more  simply.  For  this  purpose,  we  may  focus  our  attention 
on  the  vicinity  of  the  crack  tip,  since  the  effect  of  high 
stiess  concentration  is  so  dominant  there  that  the  displace¬ 
ment  field  is  almost  determined  by  the  material  strength 
alone  (Ida,  1972) .  The  known  solutions  of  dynamic  crack 
propagation  (Yoffe,  1951;  Broberg,  1961;  Kostrov,  1966) 
have  a  singularity  of  particle  v, locity  or  acceleration 
at  the  rupture  front.  This  means  that  the  initial  form 
of  source  function  corresponding  to  the  rupture  front  gives 
the  most  important  contribution  in  the  estimation  of  the 
maximum  velocitj  and  acceleration  of  the  near-field  ground 
motion.  In  this  paper,  wr  estimate  these  quantities,  based 
on  the  analysis  of  cohesive  force',  at  the  crack  tip  (Ida, 
1972)  . 


COHESIVE  FORCE  AND  PARTICLE  MOTION 

The  singularity  at  the  crack  tip  is  eliminated  by 
considerin'  an  inelastic  property  of  material  in  the  fault 
plane  near  the  tip.  According  to  Barenblatt  (1959),  the 
inelasticity  is  formulated  by  the  "cohesive  force"  that 
works  across  the  crack  tip  against  the  fracture.  Let  us 
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assume  that  the  cohesive  stress  oc  is  given  as  a  function 

of  the  displacement  discontinuity  D  across  the  fault  plane. 

The  relation  o  =  o  (D)  is  usually  called  the  "cohesive 
c  c 

force  diagram".  The  direction  of  D  relative  to  the  geo¬ 
metry  of  the  fault  piano  depends  on  whether  the  crack  is 
of  plane  shear,  tensile,  or  longitudinal-shear  type,  and 
o  (D)  is  generally  a  different  function  for  the  different 
types  of  crack.  T) o  displacement  discontinuity  D  is  an 
unknown  function  of  time  t  and  space  coordinate  x..  We 
here  choose  as  x^  the  coordinate  axis  that  is  placed 
along  the  fault  plane,  perpendicularly  to  the  crack  edge 
which  is  assumed  to  be  an  infinite  line.  The  form  of  D 
in  the  vicinity  of  the  crack  tip  is  determined  for  a  given 

a  (D) ,  as  has  been  studied  in  tie  case  of  the  longitudinal- 
c 

shear  crack  (Ida,  1972).  The  same  discussion  is  also 
applicable  to  plane  shear  and  tensile  cracks,  simply  by 
replacing  the  factor  involving  the  rupture  velocity  [see 
eqa.  (7a),  (7b)  and  (7c)]. 

To  obtain  practically  useful  expressions  for  particle 
motion,  we  introduce  a  "normalized"  cohesive  force  diagram 
0(<M,  as 


a  ^ (D)  =  ao’0(D/DQ) 


(1) 
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lhe  constants  0o  and  D0  are  the  representative  values 
of  stress  aid  displacement.  Let  us  choose  as  0o  the  stress 
of  elastic  limit  (which  is  approximately  equal  to  the 
yield  strength) ,  and  as  D0  the  displacement  required  to 
break  down  the  cohesion.  Then  we  have 

o(0)  ^  1 

3(1)  ^  0  for  t  >  1  (2) 

If  o(4>)  is  given,  a  normal  ized  displacement  discontinuity  4> 
is  determined  as  a  function  of  a  normalized  coo.dinate  X 
(Itla,  1972).  The  non-nornk*l  ized  displacement  discontinuity 
D  is  obtained  from  the  scaling  relation  [eq. (17)  of  Ida 
(1972) 1 ,  as 


D  •=  Dc4'  [-k  ( x  ^  -  ct)  J  (3) 

with  the  expression  of  X  in  terms  of  space  and  time: 

X  -  -kUj  -  ct)  (4) 

here  c  is  the  rupture  velocity  (c  >  0) ,and  k  is  defined  by 
k  =-  (327/ttu)  ( o 0/D0C )  (5) 

where  1-  is  the  rigidity,  and  Y  is  the  normalized  specific 
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surface  energy,  i.e., 

oo 

Y  =  (1/2)  |  a  ( 0 )  d4» 
o 

The  dimensionless  constant  C  in  (4)  is  a  universal 
function  of  c.  In  the  case  of  the  longitudinal-shear 
crack,  C  is  simply  (Ida,  1972): 


C 


(1  - 


cVe2)1/2 


(7a) 


where  0  is  the  shear-wave  velocity.  For  the  other  types 
of  cracks,  C  also  involves  the  longitudinal-wave  velocity 
a.  According  to  Weertman  (1969),  we  have 

C  =  4 (0/c) 2 (1  -  c2/02)_1/2[(l  -  c2/a2) (1  -  c2/02) 

-  (1  -  c2/202)  2]  (■ 

for  the  plane  shear  crack,  and 

C  =  4  (0/c)  2  (1  -  c2/a2)  1//2[(1  -  c2/a2)  (1  -  c2/02) 


(1  -  c2/202)2) 


(7c) 
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for  the  tensile  crack. 

in  fact.,  the  nature  of  cohesive  force  is  not  known 
very  well  experimentally  nor  theoretically.  Here  we 
try  to  estimate  roughly  the  effect  of  o0  an<3  Do  on  t-he 
particle  motion  wi  thout  detailed  information  on  cohesion, 
remembering  that  the  variation  of  o(<$)  is  greatly  restricted 
by  (2)  .  This  kind  of  crude  approximation  is  practically 
useful,  because  tie  parameters  o0  and  Do  can  be  estimated 
more  easily,  even  if  it  may  be  difficult  to  determine  the 
complete  form  of  cohesive  force  diagram. 

On  the  fault  plane,  the  displacement  u  is  related 
to  the  source  function  D,  as 

u  -  (1/2) D  (8) 

and  thus-  th_-  particle  velolity  u  and  acceleration  ti  are 
readily  c.itamed  frcm  the  i  unction  $ ; 

(l/2)U0i:r*<p  ' 

i)  ■=  (l/2)D0k“e^<;  " 
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where  <J>  '  and  4> "  denote  the  derivatives  of  with  respect 

to  X.  If  the  function  4>  '  (X )  has  the  maximum  value  ' 

M 

at  a  point  X  =  XM,  we  have  the  following  expression  of 
the  maximum  particle  velocity  0^: 

CiM  =  (m^'/71)  (°  0c/MC )  (9) 

Here  let  us  choose  the  orig  ine  of  X  so  that  ( 0 )  =  0. 

By  calculating  4>  (X)  for  several  models  of  cohesive  force 
diagram  a  (<{>),  it  was  shown  that  ^  '  (X)  actually  has  a 
maximum  (Ida,  1972) .  The  first  five  models  in  Table  1 
cover  the  results  in  this  calculation.  In  model  1,0' (X) 
is  infinite  at  X  =  XM,  since  a  (tj> )  is  discontinuous  at 
<J>  =  1.  In  this  model,  the  infinity  of  4>M '  is  clearly 
caused  by  the  artificial  choise  of  a(if>).  For  the  other 
models,  <f>  '  shows  a  finite  value  around  unity  in  spite 
of  the  differences  in  details  of  the  curves  a(4>).  Roughly 
speaking,  the  dimensionless  factor  1 6 Y ' /"n-  in  (9)  is 
thus  expected  to  be  of  the  order  of  the  unity.  From  (4) 
and  (5) ,  we  define  x^  and  tM  as 

xfl  =  (ttxm/32y)  (yD0/a 0C)  (10 

t w  =  X  /c  (11 

M  M 
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which  give  the  dimension  and  duration  associated  with  large 

particle  velocity.  In  the  same  manner,  we  have  the  maximum 

acceleration  ti  ,  if  <£"00  has  a  maximum  <£"  : 

M  M 

°M  *  <512Y2*M"/iT2)  (oq  c  /liDoC)2D=  (12) 

In  fact,  <£"m  is  infinite  for  all  the  models  1  through  5 
in  Table  1.  The  behavior  of  <£"(X)  is  thus  examined  in 
more  detail  in  the  next  section. 


SINGULARITY  IN  ACCELERATION 


In  the  models  1  through  5  in  Table  1,  <J>"(X)  is 

-1/2 

proportional  to  X  in  the  vicinity  of  X  =  0,  and 

not  finite  at  X  =  0.  Furthermore,  <£"(X)  is  also  in¬ 
finite  at  the  point  <£  =  i  in  some  models  (1,  2,  4  and  5). 
Of  these  two  singularities,  the  behavior  around  <£  =  1 
simply  reflects  a  more  or  less  unsmooth  connection  of  a(<J>) 
at  this  point;  we  arbitrarily  assumed  that  a(<£)  vanishes 
for  <£  >  1.  Since  the  function  <£'(X)  is  given  by  the 
Hilbert  transform  of  (X)  =  a[<£(X)]  (Ida,  1972),  <£"(X) 

is  finite  at  <£  =  1  if  a'  (<£)  is  continuous  there.  This 
condition  is  actually  satisfied  in  model  3. 


To  study  the  behavior  of  <J>(X)  in  the  vicinitv  of 
X  =  0,  we  rewrite  the  relation  between  '  (X)  and 
o32(X)  given  in  Ida  (1972),  as 

CD 

4> '  (X)  =  -(X1/2/16y)  (V  -  x)“1Y“1/2[a32(Y)  -  o32 

o 

From  (13)  ,  it  is  readily  found  that  the  following  relation 
is  necessary  so  that  <J>"(X)  may  be  finite  at  X  =  0: 

00 

Y‘3/2[a32(Y)  -  a32(0)]dy  =  0 

J  ^ 

o 

Let  us  demonstrate  through  numerical  calculation  that 
<)>"(X)  can  actually  be  finite  at  X  =  0,  if  a  suitable 
function  is  chosen  as  o  (<)>)•  We  consider  the  following 
cohesive  force  diagram,  for  e  ample; 

fa («t>  +  £  )  (<p  -  l)2  0  =  4>  =  £ 

o(4>)  = 

0  *  =  £ 

putting  a  =  2  and  £  +  i  =  1.5  in  (15),  we  keep  £ 
as  an  independent  variable.  The  examples  of  the  curve 
a(<fc)  are  given  for  several  values  of  £  in  Fig.  1. 

Since  a'(>J>)  in  (15)  is  continuous  at  <J>  =  i,  <p"  (X) 
is  finite  except  at  X  =  0  for  any  value  of  £•  The 
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distribution  of  (ji(X)  can  be  obtained  in  the  same  manner 
as  described  in  Ida  (1972).  Fig.  2  shows  the  result  for 
4>'(x).  When  £  is  larger,  the  peak  shifts  to  the  right, 
and  the  stef  p  rise  of  the  curve  around  X  =  0  is  gradually 
removed.  Finally,  at  £  =  1.37,  the  condition  (14)  is 
approximately  met,  and  if> '  (x)  increases  smmothly  at  X  =  o 
within  the  error  of  numerical  calculation.  The  distribution 
of  4>"(JC)  is  given  in  Fig.  3.  Corresponding  to  the  behavior 
of  '  (X)  ,  4>"(0)  approaches  zero  in  the  case  of  £  =  1.37. 

The  cohesive  force  diagram  for  £  =  1.37  is  denoted  by 
Model  6  in  Table  1. 

Naturally,  the  cohesive  force  diagram  that  is  physically 

realizable  must  give  finite  particle  acceleration  everywhere. 

Therefore  we  understand  that  the  condition  (14)  is  a 

physical  requirement  for  o(^).  For  such  a  (4> )  that  satisfies 

(14) ,  eq.  (12)  may  be  used  in  the  estimation  of  she  maximum 

acceleration.  The  condition  (14)  is  regarded  as  the 

criterion  to  determine  the  critical  stress  o( 0)  ,  or  o  (0), 

c 

above  which  the  deformation  is  locally  so  concentrated 
that  the  displacement  discontinuity  appears  along  a  crack 
surface  in  the  continuous  medium. 
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THE  EFFECT  OF  DISTANCE  FROM  THE  FAULT 

When  the  observation  point  is  not  situated  on  the 
fault  plane,  the  particle  motion  is  generally  affected 
by  various  factors  of  the  sources,  such  as  the  dimension 
of  the  fault  plane  and  the  starting  or  stopping  phases 
of  rupture,  other  than  the  source- time  function.  For  a 
distance  that  is  sufficiently  smaller  than  the  source 
dimension,  however,  the  stress  field  may  still  be  approx¬ 
imated  by  our  simple  2-dimensional  crack  model,  and  we  may 
estimate  how  the  particl  motion  attenuates  in  the  near 
field.  In  this  section,  we  consider  this  effect  of  dis¬ 
tance  in  the  case  of  the  longitudinal-shear  crack.  For 
the  other  types  of  cracks,  the  mathematical  treatment  will 
be  a  little  more  complicated. 

For  the  two-dimensional  longitudinal-shear  crack,  the 
displacement  u  is  given  at  an  arbitrary  point  (x^,x2) 

(Ida  and  Aki,  1972)  by 

oo 

U  =  (1/271)  (1  -  C2/B2)  1/2  x2[  [xi'2  +  (1  -  c2/32)x22]-1 

J  - 

—  00 

•D(x^  -  ct  -  x^'Jdx^'  (16) 

where  the  coordinate  x2  is  perpendicular  to  the  fault 
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plane  and  thus  |x2|  is  the  distance  from  the  fault 
plane.  Corresponding  to  (3)  and  (4),  we  introduce  the 
normalized  displacement  U  with  the  normalized  coordinate 
X2,  as 


U  =  u/Dq  (17) 

X2  =  (32Yoo/7tyDo)x2  (18) 

Because  of  the  symmetry  of  the  problem,  we  have  only  to 
consider  the  positive  side  of  X2»  The  particle  velocity 
u  and  acceleration  u  are  derived  from  U,  as 


u  =  DQkc(3U/3X) 


(19) 


m  ,2  2.^2,t/3^2. 

u  =  DQk  c  (3  U/3X  ) 


(20) 


•  II 

and  thus  the  relative  attenuations  of  u  and  u  are  given 
by  3U/3X  and  32U/DX2  as  a  function  of  X2*  For  the 
numerical  calculation,  the  following  relation  that  is 
derived  with  the  use  of  eq.  (10')  in  the  paper  of  Ida 
(1972)  is  more  convenient  than  (16): 
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au/ax  =  (1/4)  (2R) ~1/2  {(1  +  x/R)1/2 
00 

-  ( 1/ 8 y )  |  [(1  -  X/R)1/2X?  +  (1  +  X/R)1^2(Y  -  X)] 
o 

Y1/2[(Y  -  X)2  +  X22]"1o32 (Y)dY}  (21) 

where 


R 


Ui2  +  x22)1/2 


The  function  o32(Y)  denotes  the  stress  on  the  fault  plane 
(at  X2  =  0).  The  results  of  the  calculation  are  shown  in 
Figs.  4  and  5  for  model  6  of  a(<J>).  These  figures  give  the 
distribution  of  3U/3X  and  32U/3X2  as  a  function  of  X 
at  various  points  of  X2 .  The  maximum  values  of  the  nor¬ 
malized  velocity  and  acceleration  are  listed  for  various 
values  of  X2  in  Table  2.  We  may  find  that  the  acceleration 
u  decreases  rapidly  with  increasing  X2»  and  becomes  much 

smaller  than  <p  when  X_  exceeds  X  .  Therefore  the 
M2  M 

length  X M  defined  by  (10)  approximately  gives  the  distance 
within  which  we  expect  as  large  an  acceleration  as  on  the 
fault  plane,  it  is  emphasized  that  the  result  of  Table  2 
involves  only  the  effect  of  time  function  and  does  not  explain 
the  actual  ground  motion  for  very  large  distances. 
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ORIGIN  OF  STRONG  GROUND  MOTION 

It  has  been  shown  that  the  cohesive  force  diagram 
i  i)  determines  the  faulc-slip  time  function  at  the  initial 
stage,  and  yields  the  expressions  (9)  and  (12)  for  the  max¬ 
imum  velocity  uu  and  acceleration  u  on  the  fault  plane. 
For  a  rough  estimation,  we  may  drop  numerical  factors  in¬ 
volving  y,  4-  '  and  4>  " ,  and  we  have 


UM  '  (V',)C 


(22) 


(23) 


In  a  similar  manner,  we  obtain  the  following  approximate 
relations  from  (10)  and  (11): 


XM  '  (|J/00)D0 


(24) 


tM  ~  (w/oo) (D0/C) 


(25) 


Let  us  try  to  substitute  the  values  of  o  and  D  cor- 

o  o 

responding  to  the  following  two  mechanisms  of  cohesion: 


(1)  interatomic  interactions  (a 


and  D  ~  1  A)  , 
o 


and  (2)  the  gross  strength  of  rooxs  broken  by  the 


displacement  comparable  with  the  thickness  of  the  seismic 

fault  (o  ~  1  Kbar  and  D  -  10  cm)  .  The  results  of 
o  o 

the  estimation  are  given  in  Table  3.  In  the  first  case, 

we  have  very  large  velocity  and  acceleration,  but  the 

duration  or  dimension  associated  with  these  high  values 

is  too  small  to  account  for  seismic  observations.  In  the 

second  case,  we  have  the  velocity  of  1  m/s  and  the 

acceleration  of  1  g,  which  are  acceptable  values  (Drune, 

1970;  Cloud  and  Perez,  1971).  The  associated  time  scale 

t  of  0.1  second  is  also  reasonable  in  view  of  the  observed 
M 

strong  motion  spectra. 

The  results  for  uw  and  u..  in  the  second  case  agree 

M  M 

with  similar  estimates  made  by  Brune  (1970).  II-  is  emphasized, 
however,  that  the  physical  bases  for  the  estimation  are  dis- 
tinclly  different  in  the  two  studies.  In  our  treatment, 
the  maximum  particle  motion  is  determined  by  the  material 
property  concerning  the  strength  against  fracture,  while 
the  particle  velocity  or  acceleration  is  governed  by  the 
ambient  effective  tectonic  stress  in  Brune' s  model.  In 
addition,  the  predominant  frequency  of  10  cycles  is 
arbitrarily  assumed  in  Brune ' s  estimate,  but  it  was  obtained 
in  our  study.  Finally,  our  treatment  is  applicable  to  sub¬ 
sonic  rupture  propagation,  while  Brune' s  estimate  is  valid 
for  a  supersonic  propagation  which  is  not  realistic  for  a 
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spontanecus  rupture  (Ida  and  Aki ,  1972). 

The  result  in  Table  3  suggests  that  maximum  seismic 

motion  may  not  be  directly  related  to  atomic  bonding,  but 

rather  governed  by  the  bulk  strength  of  rocks.  The  values 

of  oq  and  Dq  in  the  second  case  yield  the  specific  sur- 

10  2 

face  energy  of  10  erg/cm  ,  whicu  is  much  larger  than  the 

3  2 

laboratory  measurements  10  erg/cm  for  rocks  (Brace  and 
Walsh,  1'j62).  Such  a  large  value  of  the  specific  surface 
energy  for  an  earchquake  was  also  obtained  independently 
from  the  information  on  the  time  during  which  the  rupture 
velocity  approaches  the  upper  limit  (Kikuchi  and  Takeuchi, 
1970).  This  problem  should  be  examined  more  carefully 
because  of  the  importance  of  specific  surface  energy  in 
the  fracture  phenomena. 

The  factor  C  specified  by  (6)  is  neglected  in  the 
approximate  relations  (22)  and  (23).  This  assumption  is 
valid,  unless  the  rupture  velocity  c  is  close  to  the 
critical  sound  velocity,  i . e .  ,  the  shear-wave  velocity 
for  the  longitudinal-shear  crack  and  the  Rayleigh  wave 
velocity  far  the  plane  shear  or  tensile  cracks.  In  many 
cases,  the  observed  rupture  velocity  appears  to  be  sub¬ 
stantially  smaller  than  -be  critical  value  (see  Table  3 
in  Abe,  1972),  even  for  deep-focus  earthquakes  (Fukao, 
1972).  In  such  cases,  our  rough  estimation  of  uM  and 
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(1  is  not  influenced  very  much,  even  if  the  factor  C 
M 

is  taken  into  account.  According  to  a  theoreticla  cal¬ 
culation  of  dynamic  crack  propagation,  the  rupture 
velocity  approaches  the  critical  sound  velocity  (Kostrov, 
1966) .  One  possible  explanation  of  the  discrepancy 
between  theory  and  observation  is  given  by  the  dependence 
of  Ci  on  c.  According  to  (9)  ,  u  should  be  infinitely 
large  if  c  were  the  critical  sound  velocity.  This 
suggests  that  the  energy  loss  associated  with  u,  such 
as  heat  generation , would  be  enormous  at  the  critical 
sound  velocity,  and  that  c  is  not  easily  increased 
unless  such  enormous  energy  is  supplied.  In  other  words, 
it  is  expected  that  the  theory  might  give  the  correct 
observation  of  rupture  velocity,  if  the  energy  dissipation 
is  taden  into  account. 

For  the  practical  purpose,  it  is  highly  desirable 
to  predict  the  ground  motion  in  more  detail  for  individual 
event.  At  present,  uncertain  knowledge  of  the  cohesive 
force  diagram  prevents  us  from  determining  more  accurate 
values  of  the  coefficients  in  (9)  to  (12).  Therfore 
better  estimate  cannot  be  made  only  by  giving  moic  suitable 
values  for  o0>  D0 ,  p  and  c.  An  oxpedent  way  to  overcome 
this  difficulty  might  be  to  regard  these  coefficients 
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as  empirical  parameters.  It  is  another  problem  whether 
such  large  faulting  as  earthquake  can  be  well  described 
Ijy  cohesive  force  diagram.  Further  theoretical  and 

experimental  studies  are  necessary  to  clarify  ehis  point. 
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TABLE  1 


The 

maximum  4>  '  of  the 

M 

normalized 

particle 

velocity 

4>'(X)  on  the  fault  plane 

for  given 

models  of 

cohesive 

force  diagram  a(4>) 

Model 

a  (4> )  * 

Y 

4)  ' 

XM 

1 

1 

1/2 

CD 

4.0 

2 

1-4) 

1/4 

0.683 

1.4 

3 

(1~4>)  2 

1/6 

0.909 

0.5 

4 

(1-4))  (1+24.) 

5/12 

0.693 

2.7 

S 

(l-4>2) 1/2 

tt/8 

0.712 

2.9 

6 

2  (4>  +  0.13)  (4-1. 37 )  2 

0.405 

0.581 

6.4 

*'(’*’)  0  for  4>  >  1  in  the  models  1-5,  and  for 

in  model  6. 


■!>  >1.37 
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TABLE  2 

The  maximum  values  of  the  normalized  particle  velocity 

2  2 

DU/3X  and  acceleration  3  U/3x  as  a  function  of  the 
distance  X2  from  the  fault  plane  for  model  6 


X2 

Max ( 3U/3X) 

Max(32U/3X2) 

Min(32U/3X2] 

0 

0.290 

0.16 

-0.13 

0.1 

0.266 

0.14 

-0.10 

0.5 

0.212 

0.085 

-0.048 

1.0 

0.173 

0.053 

-0.024 

5.0 

0.0873 

0.0087 

-0.0028 

10.0 

0.0624 

0.0033 

-0.0010 

50.0 

0.0280 

0.00031 

-0.000091 

100.0 


0.0198 


0.00011 


-0.000032 
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TABLE  3 

The  maximum  velocity  u„  and  acceleration  u„  for 
the  assumed  mechanisms  of  fracture 

1.  Atomic  Cohesion  2.  Strength  of  Rocks 
1  Mbar  1  Kbar 

1  A  10  cm 


fa 

o 

Assumed* 

D 

o 


Obtained 


t 


1  m/s 

1  g 

0.1  km 

0.1  s 


*  v  -  1  Mbar  and  c  ~  1  km/s  .  re  commonly  assumed. 
+  Eqs .  (22),  (23),  (24)  and  (23)  are  used. 


Fig .  1 


Fig. 


Fig. 


Fig. 


Fig . 
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FIGURE  CAPTIONS 

The  models  of  cohesive  force  diagram  given 
by  eq.  (15)  ,  in  which  a  =  2  and  i.  +  —  1  •  5 

are  assumed.  The  number  attached  to  each  curve 
denotes  the  value  of  l. 

!.  The  distribution  of  0* (X)  for  the  models  of 
cohesiva  force  in  Fig.  1.  The  number  attached 
to  each  curve  denotes  the  value  of  £. 

).  The  distribution  of  4>"(X)  for  the  model.,  of 
cohesive  force  in  Fig.  1.  The  number  attached 
to  each  curve  denotes  the  value  of  l . 

4.  The  normalized  particle  velocity  3U/aX  at 
various  distances  X2  from  the  fault  plane. 

X  may  be  regarded  as  the  coordinate  along  the 
fault,  or  the  time,  referring  to  (4). 

5.  The  normalized  particle  acceleration 

a^U/gX^  at  various  distances  X2  from  the  fault 
plane.  X  may  be  reqarded  as  the  coordinate  along 
the  fault,  or  the  time,  referring  to  (4). 


Fig. 


0.5 
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2 . 4  Generation  of  SH-Type  Motion  by  Torsion-Free  Sources 

by  S.J.  Singh 
Summary 

The  problem  of  a  spherical  cavity  in  an  infinite 
medium  is  re-examined.  The  spectral  displacement  and 
stress  fields  are  derived  when  arbitrary  tractions  are 
prescribed  over  the  surface  of  the  cavity.  This  also 
yields  the  solution  of  the  problem  of  the  release  of  pre¬ 
existing  stress  within  a  spherical  zone.  The  particular 
case  when  the  radius  of  the  cavity  is  small  in  comparison 
with  the  wavelength  under  consideration  is  discussed  in 
detail.  The  nature  of  the  tractions  are  obtained  which, 
when  applied  at  the  surface  of  the  cavity,  yield  the 

same  displacement  field  as  radiated  by  various  point 
\ 

sources.  Equivalent  body  forces  are  obtained  for  the 
release  of  pre-existing  tensile  and  shear  stresses 
within  a  spherical  zone  of  small  radius.  In  the  course 
of  the  analysis,  it  is  shown  that  an  apparently  torsion- 
free  source  is  capable  of  generating  SH-type  motion. 
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INTRODUCTION 


With  regard  to  a  given  spherical  polar  coordinate 
system,  the  motion  at  a  given  point  may  be  resolved  into  a 
toroidal  part  and  a  spheroidal  part.  For  the  toroidal 
part,  the  radial  component  of  the  displacement  and  the 
dilatation  vanish  identically.  For  the  spheroidal  part, 
the  radial  component  of  the  curl  of  the  displacement  is 
zero.  Suppose  we  shift  the  origin  of  the  coordinate  sys¬ 
tem  along  the  polar  axis.  A  motion  which  is  of  the  toroidal 
type  with  respect  to  the  old  system  will,  in  general,  contain 
both  the  toroidal  and  the  spheroidal  parts  with  reference 
to  the  new  system.  This  becomes  obvious  by  observing  that 
even  if  the  radial  component  of  the  displacement  vanishes 
in  the  old  system,  it  may  not  vanish  in  the  new  system. 
However,  if  the  motion  has  azimuthal  symmetry,  toroidal 
motion  has  only  an  azimuthal  component  while  the  azimuthal 
component  of  the  spheroidal  motion  vanishes.  In  this 
special  case,  when  the  origin  is  shifted  along  the  polar 
axis,  the  toroidal  motion  stays  toroidal  and  the  spher¬ 
oidal  motion  stays  spheroidal. 

The  realization  of  the  above  property  is  very  important 
when  considering  the  field  radiated  a  given  source.  It 
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shows  that  an  apparently  torsion-free  source  is  capable  of 
exciting  SH-type  motion. 
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SPHERICAL  CAVITY  IN  AN  INFINITE  MEDIUM 


Consider  a  homogeneous,  isotropic,  elastic  medium  of 
infinite  extent  containing  a  spherical  cavity  of  radius  a. 
Let  the  center  of  the  cavity  be  the  origin  of  a  spherical 
polar  coordinate  system  (R,0,4>).  Our  aim  is  to  calculate 
the  radiated  field  when  the  tractions  are  prescribed  over 
the  surface  of  the  cavity.  Since  the  motion  must  vanish 
at  infinity,  we  assume  the  following  expression  for  the 
displacement  at  any  point  of  the  medium 


u 


oo  a 

III 
v=c , s  1=1  m=0 


bv.nv0  (k.R) 
pm£,  mi,  3 


00(kaR)' 


(II 


where  a^,  etc.,  are  arbitrary  constants;  ka  =  w/a, 
k  =  w/6,  a )  being  the  angular  frequency  and  a,$  the  wave 

P 

velocities.  M,  N  and  L  are  the  solutions  of  the  vector 
Navicr  equation.  These  may  be  expressed  in  the  form 


^mi(kBR)  =  curl[R  hi(kBR)  YmJt(0'*)] 

=  /i  (i  +  T)£^(0, 4>)  h£(kpR)  ,  (2) 


lcBfiml<kBR)  =  CUrl  ^l<kBR> 

=  1(1+1)  £^(0,<l>)  h^tkgRj/R 

+  /m+D  ^c(©,*)(^  +  5)  h*(kBR)' 


ka^(kBR)  =  grad[h£(kaR)  ^£<0,*)] 

=  Kl{Q,*]  &h£(kaR) 

+  /m+T)  ^,(0^)  h£(kaR)/R. 


In  (2)  to  (4) ,  h  is  the  spherical  Hankel  function  of  the 

X 

second  kind  and 


Yc's(0,d>)  =  p“(cos0)(cosm$fsiwn$)  f 

mj,  x, 

=  ?R  W3'*’' 


/m+T)  3ml(e,*)  =  [e0fg  +  %STire  ^fj  Ynl(0'*)  ' 
/nr+T)  5nt(e,»)  -  (®o^Te  ^  ‘ ®' 4,)  - 


The  radial  stress  vector  is  given  by 


f  =  Xe  div  U  +  y(2-^  +  eR*curlu)  ,  (6) 


P 
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where  X.y  are  the  Lamd  parameters.  Inserting  for  u  from  (1) , 
we  find  (Ben-Menahem  and  Singh,  1968) 

oo  £ 

*R  *  E  „E,  [S^tRl/TTTiT)  2^(0.*) 

v=c,s  4=1  m=0 

+  C(R,f^,0'*)  +  ^1R)/IT!TT) 

+  Soo'W^O10'*’'  (7) 

where 


snu<R> 

*  “V  FJ.,l<rl)  °mt' 

*«t(R1 

=  2y  ( it,  (H)k.l1  ,(n) 

P  x  #  X 

<W  +  V»,j(f,Wi(,> 

WR) 

■  “fVi.2,R)  "■»  + 

2Vi,il!)  w- 

with  6  =  k  R, 
a 

n  =  kQR  and 
p 

Fl,l(lt>  ”  hl(::)  ‘  J  hUl(x)' 

,  2  (x)  t-T^2*30  "  11  h4(x)  +  lhUl(x)' 

FZ  3 (x)  =  "  !(|>21  h4<x)  +  |hUl(x)- 

'  x 


(9) 


-44- 


Let  the  tractions  at  the  surface  of  the  cavity  be  given 


by 


=  £(0,<t>)  at  R  =  a. 

R 


(10) 


We 


expand  F(0,<J>)  in  terms  of  the  vector  spherical  harmonics 


=  E  E  Ert  tav0/rrr+T)  £L(0, <t>) 

v=c,s  £=1  m=0  m£  m£ 


+  b^.(0^>  +  c^/ZTITT)  l^o,*)] 


m£  m£ 


m£ 


+  bMfSo(9'*>- 


(11) 


Given  F(0,4>),  the  expansion  coefficients  ami,  etc.  ,  can  be 
obtained  with  the  help  of  the  orthogonality  relations  for 
the  vector  spherical  harmonics  (Morse  and  Feshbach,  1953; 
p.  1900).  Equations  (7),  (10)  and  (11)  yield 


°ml<a> 


“mi' 


fVa» 


Kv  V(a) 


•  (12> 


Putting  R=a  in  (8)  and  using  (12)  ,  we  can  find  a^, 
terms  of  the  known  expansion  coefficients  am£ '  e^C ' 


etc. ,  in 
We  get 


a 


’m£ 


am»/[uk6X  Fl,l‘x)]' 


i 
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8m*  “  5sV  -  Fi,3(  =  >W'  (13) 

P  * 

YmS,  =  5ilir^  [2£(£+1)F£,l  (x)cmJl  ‘  FJt,2(x)tW]  ' 

where 

C  =  aka,  x  *  ak0,  and 
A*  =  2mi+l)FlilU)Flil(x)  -  F4,2(X)FJI,3(?)*  (14) 

This  gives  the  formal  solution  of  the  problem  of  a 
spherical  cavity  in  an  unbounded  medium  when  the  tractions 
are  prescribed  over  the  surface  of  the  cavity.  An  equivalent 
solution  was  obtained  by  Scholte  (1962).  However,  due  to  the 
difference  in  notation,  a  comparison  of  the  results  seems 
difficult. 

It  is  seen  from  (13)  that  the  original  problem  splits 

into  two  independent  problems,  one  with  bmJ^  =  c  =  0  and 

the  other  with  a  „  =  0.  The  motion  is  of  the  toroidal  type 

ms, 

in  the  former  case  and  of  the  spheroidal  type  in  the  latter 
case.  We  thus  arrive  at  the  conclusion  that  when  the  trac¬ 
tions  are  of  the  toroidal  type,  ensuing  motion  in  a  homogene¬ 
ous,  isotropic,  unbounded  medium  is  of  the  toroidal  type. 
Similarly,  spheroidal  tractions  lead  to  spheroidal  motion. 
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At  this  point,  it  may  be  noted  that  when  we  talk  about 

toroidal  or  spheroidal  type  field/  we  always  have  a  spherical 

polar  coordinate  system  in  our  minds.  In  the  above  case, 

it  is  the  (R,0,4>)  system  with  its  origin  at  the  center  of 

the  cavity.  On  the  other  hand,  the  toroidal  and  spheroidal 

motions  measured  at  the  surface  of  the  earth  ought  to  be 

with  reference  to  a  coordinate  system  having  its  origin  at 

the  center  of  the  earth.  The  coordinates  in  the  latter 

system  are  denoted  by  (r,8,<M  as  shown  in  Fig.  1.  The  two 

systems  have  the  same  polar  axis  and  rQ  is  the  distance  of 

the  center  of  the  cavity  from  the  center  of  the  earth.  The 

important  point  to  be  observed  is  that,  in  general,  a 

toroidal  or  spheroidal  type  motion  with  respect  to  the 

(R,0,<j>)  system  corresponds  to  toroidal  plus  spheroidal  type 

motion  with  respect  to  the  (r,0,4>)  system.  This  becomes 

-+•  -+• 

clear  from  the  following  transformations  which  express  M,  N 
and  L  vectors  in  the  (R,0,4>)  system  in  terms  of  these  vectors 
in  the  (r,0,4>)  system  (Wason  and  Singh,  1971)  : 


M  (knR) 

mn  p 


^mMm*(V>  lEmn(k0rO) 


U-m)kBr0^_i#i  %  <Um+l)kBr0^+1/ _ _ 

£(2£-l)  Emn  (kBr0)  "  (£+1)  (2*+3)Emn  (k0rO}  J 


£=m  J,(?+l)Emn(kBr0)  Nm£(k&r)' 


(15) 
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V<V> 


^  Nml(k6r)  tEmn(k6r0) 


( £-m)  k6rg  (A+m+l)k  r 

Emn  (k6r0)  "  7I+ITT2T+yrEmn  (k6r0}  1 


°o  k  r 

^r(TTTrEnm(k6r0)  Mm£(kBr)  ' 


l=m 


(16) 


L  ( k  R) 
mn  ot 


Z  Emn(kar0>  Lml(kar)' 


S-=TTl 


(17) 


where  r  >  rQ  and 


-*■  i 


Mml(k6r)  =  curllr  h^kgrt-^Y^U,*)], 
k6  Sml(kBr)  =  curl  "ml  (V>  ’ 


(18) 


Further,  •  Nm£^kBr^  and  Lm^kar^  are  obtained  from 

^m£,  ^kgR^  '  NmJl^k6R^  and  ^mi,^kaR^  '  resPectively ,  by  changing 

0 

R  to  r  and  0  to  0.  The  function  E  is  defined  as 

mn 


m  ,m 


(x) 


E* 

m,m+l 


(x) 


E 


i 

m,m+2 


(x) 


(2m-l)  i  I  (2U1)  {x)/xm, 

(2m+l)  !  !  (2£+l)  [(*-m)  j  (x)/xm+1 

"  jJl  +  1(x)/xm]  ,  (19) 

(2m+3)  ! I  (2£+l>  [j(£-m)  (£-m-l> j£ (x) /-m+2 
+  (m+1 )  {  j £+1  (x)  /x  -  j (x)/x  }]  , 
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where  is  the  spherical  Bessel  function  of  the  first  kind 
and 


(2m-l)!!  =  1*3*5‘  .  .  .  • (2m-l) ;  m  £  1,  (-1)!!  =  1  . 

Similar  expressions  can  be  obtained  for  other  values  of  n. 
However ,  we  will  not  need  them  here. 

From  (17)  we  note  that  an  "t  vector  in  the  (R,0,<J>) 
system  transforms  to  a  sum  of  £  vectors  in  the  (r,6,<J>) 
system.  Similarly,  if  m  =  0,  a  M  or  a  N  vector  transforms 
to  a  sum  of  M  or  N  vectors.  However,  if  m  /  0,  i.e.  if  the 
motion  is  not  symmetric  about  the  polar  axis,  a  K  or  a  N 
vector  in  the  (R,G,4>)  system  transforms  to  a  sum  of  ft  vectors 
plus  a  sum  of  N  vectors  in  t.  '  e  (r , 0  , <J> )  system.  Therefore, 
a  toroidal  or  a  spheroidal  motion  in  the  (R,0,<j>)  system  has, 
in  general,  both  the  toroidal  and  spheroidal  components  in 
the  (r,9,4>)  system. 

Once  the  field  lias  been  expressed  in  the  (r,0,<J>)  system, 
the  problem  of  satisfying  the  boundary  conditions  at  the 
surface  of  the  earth  or  any  parallel  boundary,  r  =  constant, 
becomes  straightforward.  It  can  be  easily  seen  that  if  the 
primary  field  is  of  the  toroidal  type  the  secondary  field 
will  also  be  of  the  toroidal  type.  Similarly,  for  a  primary 
field  of  the  spheroidal  type  we  get  a  secondary  field  of  the 
spheroidal  type. 
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We  may  summarize  the  above  discussion  -is  follows.  If 
the  field  due  to  a  source  with  regard  to  the  (R,0,$)  system 
is  axially  symmetric  and  is  of  the  toroidal  or  spheroidal 
type  then  the  total  field  due  to  such  a  source  in  a  radially 
heterogeneous  earth  will  be  of  the  toroidal  or  spheroidal 
type,  respectively,  even  when  referred  to  the  (r,0,4>)  system. 
But  if  the  field  due  to  the  source  is  not  axially  symmetric 
and  is  of  the  toroidal  or  spheroidal  kind  with  regard  to  the 
(R, 0 , 4>)  system,  the  total  field  when  the  source  is  placed 
in  a  radially  heterogeneous  earth  has  both  the  toroidal  and 
spheroidal  components  with  regard  to  the  (r , 0  , <*> )  system. 


THE  CASE  OF  SMALL  CAVITY 


We  next  assume  that  the  radius  of  the  cavity  is  small 
in  comparison  with  the  wavelength  under  consideration  so 
that  C  <<  1,  X  <<:  1.  We  can  then  use  the  following  expansion 
of  the  spherical  Hankel  function  about  the  origin: 


Vx) 


1  ,._A  1  .M2 

(2U1)  !  1  2  (2A+3)X 


+  TTTTTI+iy  (2£+5Tx£+4  +  *"1 

-2.-1  1  -2.+  1 
+  i(2A-l)!![x  *  +  2(2Xriy* 


1  -£+3  ,  n 

+  2.4  in-\)  Ul-T)x  +  •  •  • J  • 


(20) 


We  shall  be  considering  only  such  tractions  as  could  be 
expressed  in  terms  of  the  spherical  harmonics  of  degree 
i  <  2.  Using  (9),  (13),  (14)  and  (20),  we  obtain  the  follow¬ 
ing  limiting  values  for  a etc.  for  M  2  as  £  and  x  tend 
to  zero: 


II 

o 

o 

1  J)  —  w 

~T y  a  ka  b00; 

(21) 

°m,  1 

i  3,  2 

3y  a  k6  am,  1 ' 

Bm,l 

(22) 

ii 

i"H 

(B/o.)3  Bmili 
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a 


’m ,  2 


m,2 


i  4.  3 

12y  B  am,2’ 


l 

V 


X  +  2y  3.  2,,  ^  _ 

9X  +  14 y J  a  kfl  ^b™  “J  +  3c~  ^  ' 


B  m,2  m,2 


(23) 


V2  '  2<8/“>  «m,2 


Inserting  the  above  values  in  (1)  and  writing 


S  -  i  l  z  i\  +  s= 

«  0_1  _n  00' 

v=c,s  &=1  m=0 


(24) 


we  find 


■+c 

U00 


-4*  aVbS0  200(ltc,R>' 


(25) 


a3k  2aV  ,  Mv  (k„R) 
m,l  3y  B  m,l  m,l  B 


+  a2k„  (bv  .  +  2cv  .)  [Nv  .  (k„R) 
3y  B  m,l  m,l  m,l  B 


+  (B/0)3  ^,l(kaR)1’ 


(26) 


-*-v 

U  o 
m ,  2 


yi-  a4k  3aV  7  Mv  . (kfiR) 

12y  B  m,2  m,2  6 


i  X  +  2y  _3.  2  ,,  v  ,  v  . 

y  9X  +  14y  a  kB  vbm,2  3cm,2) 


m 


IN. 


4  +V 


,  2  (kgR)  +  2(B/a)n  2(kaR)]  •  (27) 


Equation  (24)  gives  the  displacement  field  in  an 
unbounded  medium  containing  a  small  spherical  cavity,  of 
radius  a,  the  boundary  of  which  has  prescribed  surface  trac¬ 
tions  F(0,  $).  Equations  (25)  to  (27)  are  similar  to  the 
expressions  for  the  displacement  field  due  to  various  point 
sources  situated  in  an  infinite  medium.  We  shall  derive  the 
relationship  between  the  seismic  moment  of  a  point  source 
and  the  tractions  at  the  surface  of  the  cavity  so  that  the 
two  yield  identical  radiations.  It  has  been  shown  by  Singh 
—  •  (1972)  that  the  displacement  field  due  to  a  point 

source  in  a  homogeneous,  isotropic,  infinite  medium  may  be 
expressed  in  terms  of  the  vectors  M,  N  and  L.  In  the  follow¬ 
ing,  the  results  for  point  sources  have  been  taken  from 
Table  3  of  Singh  et  al.  (1972) . 


(i)  Center  of  compression 

A  center  of  compression  is  equivalent  to  three  equal 
mutually  perpendicular  dipoles.  The  seismic  moment  of  a 
dipole  is  defined  as  the  product  of  the  force  and  arm-length. 
The  seismic  moment  of  a  center  of  compression  may  be  defined 
as  the  seismic  moment  of  any  one  of  the  component  dipoles 
(Aki  and  Tsai, 1972).  Let  the  seismic  moment  of  the  center 
of  compression  be  Mg.  We  then  have 


u  = 


lMokg 

4iry 


(*> 

a 


<V> 


(28) 
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Comparing  (25)  and  (28) ,  we  find 


00 


Mo 

tt  a 


y 

X  +  2y 


(29) 


It  is  obvious  from  (11)  and  (29)  that,  in  this  case,  the 
surface  traction  is  simply  a  pressure.  The  relation  between 
the  magnitude  of  the  pressure,  pQ,  and  the  siesmic  moment, 

M  ,  of  the  equivalent  center  of  compression  is 


Mq  =  TraJ  pQ  (X  +  2P)/y . 


(30) 


This  relation  has  been  given  earlier  by  Aki  and  Tsai  (1972) . 


(ii)  Center  of  rotation 

If  the  rotation  is  about  the  z-axis,  we  have 


-*■  iMo’K8  •>c 

u  =  --TO—M01(kBR) 


(31) 


Comparing  with  (26) ,  we  get 


3M 

4TTa3 


(32) 


the  rest  of  the  coefficients  being  zero. 

Equation  (11)  shows  that,  for  a  center  of  rotation  about 
the  z-axis,  of  seismic  moment  , 
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Trr  =  0,  tdR  =  0*  T0.  =  -^sin6» 


RS 


R  (p 


O' 


where 


M0  =  T  lr,1^T0‘ 


(iii)  Single  couple 

Tor  a  couple  of  moment  in  the  xy-plane  with  its 
forces  parallel  to  the  x-a*is, 


> 

u 


f  .  jts  oAVs  , 

45rTi"i"6ll0i  +  N22  +  2^a)  L22^* 


(33) 


Comparing  (26)  ,  (27)  and  (3  3)  ,  we  obtain 
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3  r  +  2y  J 


(34) 


If  the  for  ;o;>  of  the  couple  are  parallel  to  the  y-axis,  we 
have  instead 


“  ’  -  Wu-[6fiSl  +  «S22  +  2(§,4J221- 


(35) 
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From  (34)  ,  we  get  b^2  +  ^c22'  not  b 22  anc^  c22 

separately.  This  implies  that  the  solution  is  not  unique. 
Corresponding  to  a  given  source,  we  can  find  more  than  one 
distribution  of  tractions  over  the  surface  of  a  small 
cavity  which  would  yield  the  same  radiations  as  the  given 
source.  This  lack  of  uniqueness  is  obviously  caused  by 
the  limiting  process.  For  a  cavity  of  radius  which  may 
not  be  small,  we  have  seen  in  (7),  (8)  and  (13)  that  a 
given  displacement  field  gives  a  definite  stress  field  on 
the  surface  of  the  cavity  and,  conversely,  corresponding  to 
a  given  distribution  of  stresses  at  the  surface  of  the 
cavity,  there  is  a  definite  unique  displacement  field. 


(iv)  Double  couple 

Defining  a  double  couple  of  seismic  moment  Mq  in  the 
xy-plane  as  the  sum  of  the  two  single  couples  whose  dis¬ 
placement  fields  are  given  by  (33)  and  (35)  ,  we  have 


u  = 


and 


bh  +  3c22 
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(36) 
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(v)  Dipole 

For  a  dipole  along  the  z-axis  of  seismic  moment  Mg 


u  =  - 


iMoV 

l2rrvi 


'^02 


-  (tSo  -  2l02n-  (38) 


As  before,  we  find 


X  +  2yJ  ' 


M, 


12fra' 


91  +  14y 
X  +  2y  j 


(39) 


However,  if  the  dipole  is  along  the  x-axis,  we  have 


-*•  lMoke  roac  $c 
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(40) 


For  this  case. 
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We  have  obtained  above  tractions  or  a  set  of  tractions 
which  when  applied  over  the  surface  of  a  small  cavity  in  an 
infinite  medium  give  the  same  field  as  radiated  by  a  point 
source  in  an  infinite  medium  in  the  absence  of  the  cavity. 


RELEASE  OF  STRESS  IN  A  SPHERICAL  ZONE 


Consider  a  homogeneous,  isotropic,  elastic,  prestressed 
medium.  Suppose  that  the  stress  is  released  within  a  spherical 
zone  of  radius  a.  Our  aim  is  to  calculate  the  additional 
displacement  field  caused  by  the  release  of  the  pre-existing 
stress.  Denote  the  initial  radial  stress  vector  at  any 
point  by  TR^  and  let 

Tr°  =  -  F  (0 ,4> )  at  R  =  a.  (42) 


If  the  additional  displacement  field  is  expressed  as  in  (1) 
and  the  corresponding  stress  vector  as  in  (7),  then  we  must 
have 

■fR  +  $R°  =  0  at  R  =  a.  (43) 
Equations  (42)  and  (43)  yield 


$_  =  -  ?_°  =  ?(©,$)  at  R  =  a.  (44) 

K  K 

Comparing  (10)  and  (44)  we  note  that  the  additional  field 
due  to  the  release  of  the  stress  in  a  spherical  zone  is 
identical  with  the  field  due  to  a  spherical  cavity  in  an 
infinite  medium  with  prescribed  surface  tractions. 


i 
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Consequently ,  we  can  use  here  the  results  derived  before 
for  a  spherical  cavity,  including  those  for  a  cavity  of 
small  radius. 

Consider  the  case  in  which  the  initial  stress  is  a 
simple  tension  of  magnitude  Tq  parallel  to  the  axis  of  x. 
We  have 


T°  =  T  ,  T°  =  T°  =  T0  =  T°  =  T°  =  0, 
xx  O'  yy  zz  yz  zx  xy 


(45) 


Equations  (42)  and  (45)  yield 
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using  (5).  Comparing  (11)  and  (46),  we  get 
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The  coefficients 

not  mentioned  are  zero.  Inserting 

(13),  one  finds  the  coefficients  am^,  etc. 
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(48) 


Inserting  into  (1)  •  we  obtain  the  displacement  field 
due  to  the  stress  release 
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„c  ,-tc  l=fcC  .  ,  c 

"  sG2(N02  "  z  22  +  Y00iJ00 
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(49) 


Once  again,  we  consider  the  case  when  the 
soherical  zone,  a,  is  small  in  comparison 
length  of  interest.  Therefore,  C  «  1,  X 
(14),  (20,  and  (48)  yield 


radius  of  the 
with  the  wave- 
<<  1 ,  and  (9) , 
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Equation  (49)  now  becomes 
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Comparing  the  last  equation  with  (28)  and  (40),  we  note 
that  the  displacement  field  (51)  is  identical  with  the  field 
due  to  a  dipole  in  the  x-direction  of  seismic  moment 


M, 


=  20ira  t, 
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9T 


2y  ) 
14yJ 


(52) 


plus  a  center  of  compression  of  seismic  moment 
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(X  +  2y)  (3X  -  2y ) 
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(53) 


We  thus  see  that  the  displacement  field  due  to  the 
release  of  a  pre-existing  stress,  in  the  form  of  simple 
tension  parallel  to  the  x-axis,  within  a  small  spherical 
zone,  is  identical  with  the  displacement  field  due  to  a 
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tension  dipole  in  the  x-direction  and  a  center  of  compres¬ 
sion.  Similarly,  if  the  pre-existing  stress  is  a  simple 
pressure  in  the  y-direction,  the  equivalent  force  system 
is  a  pressure  dipole  in  the  y-direction  plus  a  center  of 
dilatation.  Combining  the  two  results,  we  _ind  that  if  the 
initial  stress  is  a  tension  Tq  in  the  x-direction  plus  a 
pressure  Tq  in  the  y-direction,  the  equivalent  force  system 
^s  a  tension  dipole  in  the  x-direction  together  with  a 
pressure  dipole  in  the  y-direction,  each  of  seismic  moment 
Mq  given  by  (52) . 

Suppose  a  new  coordinate  system  [x*  ,y',z')  is  obtained 
by  rotating  the  (x,y,z)  system  about  the  z-axis  through  an 
angle  -r/4 .  In  that  case 


(e  -  e  )/n, 
x  y  '  ' 

(5X  +  £y>//*. 

This  yields 


->  •+  -+•  ->• 
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(54) 


(55) 


Tne  last  result  shows  that  a  tension  Tq  in  the  x-direction 
and  a  pressure  t  ^  in  the  y-direction  are  equivalent  to  a 
shear  with  reference  to  the  x'y'  directions.  Further, 
we  know  tnat  a  tension  dipole  Mq  and  a  pressure  dipole  Mq 
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acting  along  mutually  perpendicular  directions  are  together 
equivalent  to  a  double  couple  Mq.  We,  therefore,  conclude 
that  the  displacement  field  due  to  the  release  of  pre¬ 
existing  shearing  stress  Tq  within  a  small  spherical  zone 
of  radius  a  is  identical  with  the  displacement  field  due 
to  a  double  couple  of  seismic  moment  Mq ,  where 


M  on  ,3  1  +  A 

Mq  207ra  Tq  9^  + 


This  result  could  have  been  derived  directly  by  assuming  an 
initial  stress  of  the  form 


0  °  =  0  =  0  =  0=t°=0 

Txy  "  T0'  T xx  Tyy  Tzz  Tyz  Tzx 


(57) 


For  the  Poisson  case  (X  =  y) ,  (56)  becomes 


60  3 

U  ta  Tq- 


(58) 


Equation  (58)  has  been  obtained  earlier  by  Aki  et  al .  (1969) 

and  by  Aki  and  Tsa i  (1972),  beginning  with  the  results  of 
Honda  (1960,  1962).  Since  Honda  gives  only  the  far-field 
displacements,  Aki  et:  al .  (1969)  and  Aki  and  Thai  (1972) 

state  that  relation  (53)  holds  only  for  the  far-field.  How¬ 
ever,  as  we  nave  seen  above,  the  relationship  applies  for 
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the  totax  exact-field.  The  only  assumption  is  that  the 
radius  of  the  spnerical  zone  is  small,  in  comparison  with 
the  wavelength  under  consideration. 
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2 . 5  On  the  Disturbance  Due  to  a  Spherical  Distortional  Pulse 
in  an  Elastic  Medium  by  S.J.  Singh  and  M.  Rosenman* 
(Abstract) 

Theoretical  expressions  are  derived  for  the  displace¬ 
ment,  velocity  and  stress  in  the  time  domain  induced  by  an 
axially  symmetric  shearing  stress  applied  at  the  inner  sur¬ 
face  of  a  spherical  cavity  in  a  homogeneous,  isotropic, 
elastic  medium  of  infinite  extent.  Theoretical  seismograms 
are  computed  for  a  step  source  and  for  three  sources  with 
exponential  decay  in  time.  A  satisfactory  time-dependence 
of  the  source  can  be  obtained  by  combining  the  step  source 
with  one  or  more  exponentially  decaying  sources. 

Address;  Department  of  Geology,  Harvard  University, 
Cambridge,  Massachusetts. 

2 • 6  A_^£herical  Cavity  in  a  Micropolar  Elastic  Medium  and 
Related  Problems ,  by  S.J.  Singh  (Abstract) 

The  problem  of  a  spherical  cavity  in  an  infinite, 
linear,  isotropic,  micropolar  elastic  medium  is  considered. 
The  spectral  displacement  and  stress  fields  are  obtained 
when  arbitrary  tractions  and  couples  are  prescribed  over 
the  surface  of  the  cavity.  It  is  found  that,  as  in  the 
elastic  case,  the  original  problem  splits  into  two 
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independent  problems  corresponding  to  spheroidal  and  tor¬ 
oidal  motions,  respectively.  In  the  case  of  purely  radial 
static  or  dynamic  displacements,  there  is  no  microrotation. 
The  solutions  of  such  micropolar  elastic  problems  can  be 
obtained  from  the  corresponding  elastic  solutions  on  re¬ 
placing  y  by  y  +  ^-k  .  This  correspondence  principle  is 
used  to  derive  the  micropolar-elastic  solutions  of  several 
problems  involving  radial  displacements  only. 

2 . 7  Quasi-Static  Deformation  of  a  Viscoelastic  Half-Space 

by  a  Displacement  Dislocation  by  S.H.  Singh  and 

M.  Rosenman  (Abstract) 

The  correspondence  principle  is  used  to  get  the  Laplace 
transformed  solution  of  the  problem  of  the  quasi-static 
deformation  of  a  viscoelastic  half-space  by  a  shear  dis¬ 
placement  dislocation  from  the  corresponding  elastic  results. 
The  transformed  solution  is  inverted  for  the  Voigt  and  the 
Maxwell  viscoelastic  models.  It  is  shown  that,  for  a  verti¬ 
cal  dip-slip  fault,  the  surface  displacements  for  the  visco¬ 
elastic  case  are  identical  with  the  elastic  displacements. 

In  the  case  of  a  vertical  strike-slip  fault,  detailed 
numerical  results  are  obtained  for  both  a  point  source  and 
a  finite  rectangular  source.  It  is  found  that  the  results 
for  the  viscoelastic  models  differ  significantly  from  the 
corresponding  elastic  results. 
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2 • 8  Quasi-Static  Strains  and  Tilts  Due  to  Faulting  in  a 

Viscoelastic  Half-space  by  M.  Rosenman  and  S.J.  Singh 
(Abstract) 

Recently  derived  quasi-static  surface  displacements 
resulting  from  a  finite  rectangular  vertical  strike-slip 
fault  in  a  viscoelastic  half-space  are  used  to  derive 
the  surface  strains  and  tilts  for  both  Voigt  and  Maxwell 
viscoelastic  models.  Contour  maps  are  obtained  for  various 
strain  and  tilt  components.  The  variation  with  time  and 
epicentral  distance  is  studied  in  some  representative  cases. 
Detailed  numerical  calculations  reveal  significant  differences 
between  the  viscoelastic  and  the  elastic  results  in  the  case 
of  a  vertical  strike-slip  fault.  This  contrasts  with  the 
results  for  a  vertical  dip-slip  fault  in  a  uniform  half¬ 
space,  where  the  surface  strains  and  tilts  for  the  visco¬ 
elastic  and  the  elastic  models  are  identical. 
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3.  EARTH  STRUCTURE  AND  PATH  EFFECTS 

3,1  Regionalized  Earth  Models  from  Linear  Programming 

Methods  by  C.E.  Johnson  (Abstract) 

This  study  is  concerned  with  the  development  of  pos¬ 
sible  models  of  the  internal  structure  of  the  earth  con¬ 
sistent  with  a  given  set  of  observed  data.  A  two-stage 
linear  programming  procedure  was  used  together  with  an 
assumed  parameterization  to  obtain  an  explicit  envelope 
of  possible  shear  velocity  and  density  values  in  the  mantle 
and  core.  This  envelope  is  determined  separately  for  oceanic 
shield,  and  tectonic  regions  of  the  upper  mantle.  The  data 
used  in  this  study  consist  of  the  mass  and  moment  of  inertia 
of  the  earth,  periods  of  free  oscillations,  including  re¬ 
cently  available  overtones,  regionalized  phase  and  group 
velocities  of  Rayleigh  waves,  and  phase  velocities  of  Love 
waves.  The  results  constrain  the  variations  of  density  and 
shear  velocity  in  the  lower  mantle  to  within  about  1.5% 
from  the  center  of  the  envelope.  The  density  just  below 
the  mantle-core  boundary  was  found  to  lie  between  9.79  and 
9.86  grams/cc.  A  rigid  core  was  needed  to  satisfy  the  over¬ 
tone  data  with  a  shear  velocity  between  3.35  and  3.52  km/sec. 
The  radius  of  the  mantle-core  boundary  was  found  to  lie  be¬ 
tween  3476.38  and  3486.42  kilometers. 


Excellent  agreement 
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with  recent  travel  time  studies  of  body  waves  was  found  for 
shear  velocity  in  the  lower  mantle  and  for  the  radius  of 
the  mantle-core  boundary.  Geophysical  and  petrological 
interpretations  based  on  these  results  are  discussed. 

3.2  Range  of  Earth  Structure  Nonuniqueness  Implied 

by  Body  Wave  Observations  by  R.A.  Wiggins,  G.A.  McMechan 
and  M.N.  Toksdz  (Abstract) 

The  Herglotz-Wiechert  integral  for  the  direct 
inversion  of  ray  parameter  versus  distance  curves  can  be 
manipulated  to  find  the  envelope  of  all  possible  models 
consistent  with  geometrical  body  wave  observations 
(travel  time  and  ray  parameter  versus  distance) .  Such 
an  extremal  inversion  approach  has  been  used  to  find 
the  uncertainty  bounds  for  the  velocity  structure  in 
the  mantle  and  core.  We  find,  for  example,  that  there 
is  an  uncertainty  of  ±40  km  in  the  radius  of  the  inner 
core  boundary,  ±18  km  at  the  core-mantle  boundary,  and 
±35  km  at  the  435-km  transition  zone.  The  velocity 
uncertainty  is  about  ±0.08  km/sec  for  P  and  S  waves  in 
the  lower  mantle  and  about  ±0.20  km/sec  in  the  core. 
Experiments  with  various  combinations  of  ray  tapes  in 
the  core  indicate  that  rather  crude  observations  of 
SKKS-SKS  travel  times  confine  the  range  of  possible 
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models  far  more  dramatically  than  do  the  most  precise 
estimates  of  PmKP  travel  times.  Comparisons  of  results 
from  extremal  inversion  and  linearized  perturbation 
inversions  indicate  that  body  wave  behavior  is  too 
strongly  nonlinear  for  linearized  schemes  to  be  effective 
for  predicting  uncertainty. 

3 . 3  A  Travel-Time  Study  Using  Deep-Focus  Earthquakes 

by  M.K.  Sengupta  (Abstract) 

Revision  of  the  Jef f reys-Bullen  (J-B)  tables  has  been 
made  from  the  analysis  of  travel-times  of  deep  earthquakes. 
Deep  events  (450  to  600  km)  were  used  in  order  to  avoid 
the  source  bias  caused  by  a  downgoing  slab.  The  absolute 
values  of  travel-times  have  been  determined  from  the 
Nevada  Test  Site  explosion  data  for  which  the  upper  mantle 
velocity  structure  near  the  source  is  known  and  could  be 
corrected  for . 

In  this  analysis,  station  errors  and  the  systematic 
error  of  the  J-B  tables  are,  in  general,  similar  to 
other  works.  However,  a  smaller  scattering  of  the  data 
suggests  that  the  results  of  this  study  may  be  more 
reliable.  Prediction  of  travel-times  from  deep  events, 
barring  a  d.c.  term,  could  be  made  from  our  times  and 
station  corrections  with  an  error  comparable  to  the  reading 


error. 
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Residual  sphere  plots  of  data  suggest  that  anomalous, 
high  velocity  structures  of  the  mantle  beneath  island 
arcs  might  exist  even  beyond  the  region  of  deep  earth¬ 
quakes  (e.g.,  in  the  Solomon  islands).  Also,  there  is  a 
suggestion  of  lateral  heterogeneity  at  the  deep  mantle  as 
evidenced  from  these  plots  and  the  anomalous  travel-time 
variations  between  different  source  zones  beyond  80° 
of  distance.  Stations  in  western  North  America  were 
found  to  show  different  residuals  for  different  source 
regions  in  different  azimuths  reflecting  a  complicated 
velocity  structure  underneath  the  stations. 

3 . 4  Seismic  Travel-Time  Evidence  for  Lateral  Inhomogeneity 
in  the  Deep  Mantle  by  B.R.  Julian  and  M.K.  Sengupta 
(Abstract) 

Regional  variations  in  travel  times  at  distances 
greater  than  about  70°  have  been  found  in  a  study  of 
P  waves  from  deep  focus  earthquakes.  These  data  can  be 
explained  satisfactorily  only  in  terms  of  large-scale 
lateral  heterogeneity  in  the  lowest  few  hundred  km  of 
the  mantle,  this  region  being  more  heterogenous  than 
that  which  lies  above  it.  The  travel  time  varies  by 
more  than  a  second,  indicating  at  least  a  1%  variation 
in  the  P  wave  velocity  in  the  deep  mantle.  Among  other 
results,  the  data  indicate  a  pronounced  lateral  variation 
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beneath  the  Hawaiian  islands,  the  velocities  being  high 
to  the  northwest  of  Hawaii  and  low  in  the  vicinity  of 
the  islands.  It  also  appears  that  low  velocities  may 
be  characteristic  of  island  arcs  at  depths  greater  than 
1000  km. 

3.5  Scattering  of  P  Waves  Under  the  Montana  LASA 

by  K.  Aki  (Abstract) 

The  variations  of  amplitude  and  time-delay  of  teJe- 
seismic  P  waves  across  the  Montana  LASA  were  interpreted 
as  due  to  scattering  by  a  random  inhomogeneity  in  the 
earth's  crust  under  the  array.  The  prediction  of  the 
Chernov  theory  explains  well  the  observed  statistical 
properties  of  P  waves  with  frequency  0.5  cps .  The  in¬ 
homogeneity  under  LASA  has  a  correlation  distance  of  about 
10  km,  with  a  fractional  RMS  velocity  fluctuation  of 
4%,  extending  to  a  depth  of  about  60  km.  The  turbidity 
coefficient  under  LASA  at  0.5  cps  is  0.008  km-1,  which 
is  much  greater  than  the  values  (10_^  ~  10~4km_1  at  5  ~  10 
cps)  observed  for  refracted  waves  in  the  crust  and  upper 
mantle  in  the  U.S.S.R.  by  Nikolayev  and  his  colleagues. 

The  scattering  under  LASA  is  so  strong  that  the  condition 
for  the  Born  ipproximation  is  violated  for  frequencies 
higher  than  0.5  cps.  Accordingly,  the  observed  statistical 
properties  of  1  cps  waves  show  systematic  departure  from 
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the  prediction  of  the  Chernov  theory. 

3 . 6  Shear-Wave  Attenuation  and  Melting  Beneath  the 

Mid-Atlantic  Ridge  by  S.C.  Solomon 
Summary 

Because  the  attenuation  of  seismic  waves  is  sensitive 
to  variations  in  temperature  and  to  partial  melting,  the 
mapping  of  seismic  Q  beneath  the  mid-ocean  ridge  systems 
is  a  useful  tool  to  outline  boundaries  between  lithosphere 
and  asthenosphere  and  to  constrain  the  mechanics  of  the 
intrusion  process.  Our  approach  is  to  measure  the  differ¬ 
ential  attenuation  of  long-period  shear  waves,  using  a 
spectral  ratio  technique,  from  earthquakes  on  the  ridge 
and  to  look  for  variations  in  attenuation  with  propagation 
direction.  We  correct  for  proprgation  distance  and,  where 
known,  the  upper  mantle  attenuation  beneath  the  receiving 
stations . 

The  azimuthal  dependence  of  attenuation  of  S  waves 
from  an  earthquake  offset  from  the  ridge  axis  on  a  trans¬ 
form  fault  indicates  the  existence  of  a  low-Q  zone,  no 
wider  than  100  km  and  confined  to  depths  shallower  than 
about  100  to  150  km,  beneath  the  crest  of  the  mid-Atlantic 
ridge.  The  absence  of  appreciable  azimuthal  variation  in 
shear-wave  attenuation  for  an  earthquake  on  the  ridge 
crest  suggests  that  the  low-Q  zone  is  at  least  50  km  wide. 
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Q  within  such  a  zone  must  be  10  or  less  for  long-period 
S  waves.  The  most  likely  explanation  of  such  a  low-Q 
zone  of  limited  spatial  extent  and  with  sharply  defined 
boundaries  is  that  the  zone  is  a  region  of  extensive 
partial  melting,  probably  at  temperatures  in  excess  of 
the  anhydrous  solidus  of  mantle  material.  Such  a  region 
of  large  melt  concentration  is  consistent  with  the 
chemistry  of  rocks  from  the  mid-ocean  ridges  and  with 
models  of  the  temperature  field  derived  from  numerical 
calculations  of  flow  beneath  spreading  centers. 
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INTRODUCTION 

The  upper  mantle  beneath  mid— ocean  ridges  and  other 
spreading  centers  is  anomalously  hot,  most  likely  because 
the  ascending  limb  of  some  form  of  convection  cell  is 
centered  more  or  less  beneath  the  ridge  axis.  Portions  of 
the  mantle  beneath  a  spreading  ridge  must  further  be  par¬ 
tially  n j1 ten ,  if  for  no  other  reason  than  to  segregate 
the  basaltic  magmas  needed  to  form  the  oceanic  crust.  High 
temperatures  and  partial  melting  can  have  a  profound  effect 
on  the  velocity  of  seismic-wave  propagation  and  especially 
on  the  seismic  attenuation  or  the  quality  factor  Q.  Beneath 
mid-ocean  ridges,  a  zone  of  unusually  low  velocity  and  low 
Q  extends  approximately  from  the  base  of  the  crust  [Le  Pichon 
et  al . ,  1965;  Molnar  and  Oliver,  1969;  Keen  and  Tramontini , 
1970],  at  least  in  the  immediate  vicinity  of  the  ridge  crest, 
to  several  hundred  kilometers  depth  [Francis ,  1969;  Knopof f 
et  al.  ,  1970;  Weidner,  1972;  Forsyth,  1972;  Taylor,  1972]. 
This  zone  is  coi.-i.’only  identified  with  a  partially  molten 
asthenosphe re ,  whose  accentuated  features  and  shallow  top 
are  a  consequence  of  thi  elevation  of  isotherms  by  convective 
upwel  1  i  n<i . 

Rotter  quantitative  information  on  the  spatial  variation 
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of  (J  beneath  spreading  mid— ocean  ridges  can  therefore  serve 
as  a  valuable  constraint  on  the  mantle  temperature  field  and 
the  process  of  lithosphere  generation.  While  the  most  de¬ 
tailed  such  information  will  probably  result  from  short- 
range  refraction  experiments  using  ocean-bottom  senseis, 
we  are  limited  at  present  to  land-based  observations.  The 
approach  in  this  paper  will  be  to  measure  the  amplitude 
spectra  of  long-period,  teleseismic  shear  waves  generated  by 
earthquakes  on  a  mid-ocean  ridge,  to  look  for  changes  in  the 
differential  attenuation  of  the  waves  with  propagation  direc¬ 
tion,  and  to  relate  such  changes  to  spatial  variations  of  Q 
beneath  the  ridge.  The  most  important  result  of  this  study 
is  the  striking  evidence  for  a  narrow,  shallow  region  of 
very  low  Q  lying  beneath  the  crest  of  the  mid-Atlantic  ridge. 
This  low-Q  zone  is  clearly  associated  with  melting,  and 
probably  represents  the  region  in  whi.ch  temperatures  exceed 
the  dry  solidus  of  mantle  material. 


SOURCE  OF  DATA 

Wo  shall  present  in  this  paper  the  differential  attenu- 
at  ion  ot  waves  from  two  earthquakes  on  the  mid-Atlantic 


i  i  dge  . 


All  data  foi  fins  work  are  taken  from  recordings  on 
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horizontal,  long-period  seismometers  of  the  World-Wide 
Standard  Seismograph  Network.  The  mid-ocean  ridge  in  the 
north  Atlantic  is  an  ideal  spreading  center  to  be  studied 
using  toleseismic  body  waves  because  of  the  excellent  azi¬ 
muthal  coverage  of  standard  stations.  Pertinent  information 
on  the  two  earthquakes  treated  is  given  in  Table  1.  Their 
epicenters  are  shown  superimposed  on  a  seismicity  map  for 
the  north  Atlantic  in  Fiqure  1. 

The  earthquake  of  2  June  1965  was  located  on  the  crest 
of  the  mid-Atlantic  ridge.  The  focal  mechanism  is  one  of 
normal  faulting,  with  the  inferred  tensional  pre-stress 
axis  approximately  horizontal  and  perpendicular  to  the 
strike  of  the  ridge  [Sykes ,  1970a]  .  We  shall  contrast  below 
the  attenuation  of  shear  waves  that  have  propagated  more 
nearly  parallel  to  the  ridge  axis  with  those  that  propa¬ 
gated  roughly  perpendicular  to  the  axis.  Of  some  note  is 
the  eastward  offset  in  the  ridge  crest  immediately  to  the 
south  of  the  epicenter.  Further  mention  will  be  made  of 
this  point  later. 

The  earthquake  of  13  February  1967  was  located  on  the 
Gibbs  fracture  zone,  a  prominent  transform  fault  that  offsets 
two  port  ions  of  the  mid-Atlantic  ridge  crest  by  over  350  km 
(Fleming  et  ,d_.  ,  1  970].  The  focal  mechanism  derived  from 

I'-wave  first  motions  is  shown  in  Figure  2;  clearly  right- 
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lateral  strike-slip  motion  on  a  nearly  vertical  fault  is 
indicated.  The  strike  of  the  east-west  nodal  plane  is 
identical  to  the  strike  of  the  fracture  zone  [Fleming  et  al . , 
1970].  The  earthquake  epicenter  is  roughly  60  km  east  of 
the  northern  ridge  crest  segment.  We  shall  pay  particular 
attention  below  to  the  contrast  in  attenuation  between  shear 
waves  that  have  propagated  beneath  the  ridge  crest  and  those 
that  have  traveled  beneath  the  older  ocean  floor  to  the  south. 
We  might  remark  that  the  Gibbs  fracture  zone  is  somewhat  unu¬ 
sual.  It  is  actually  a  double  fracture  zone  and  apparently 
marks  a  boundary  between  regions  with  slightly  different 
spreading  directions  [Fleming  et  al . ,  1970] .  Added  to  the 
confusion  is  that  spreading  directions  predicted  for  this 
region  of  the  north  Atlantic  from  Europe-North  American 
spreading  poles  [Le  Pichon,  1968;  Chase ,  1972]  do  not  match 
the  trend  of  the  fracture  zone,  though  these  two  determinations 
of  spreading  direction  differ  by  17  degrees  and  bracket  the 
strike  direction  of  the  transform  fault. 
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DIFFERENTIAL  ATTENUATION 

Background 

A  seismic  signal  contains  information  on  both  the 
source  of  the  signal  and  the  medium  through  which  it  has 
propagated .  To  isolate  the  effect  of  seismic  attenuation, 
we  appeal  to  the  technique  of  body-wave  equalization  [Ben- 
Monahem  — -  a -*  *  1965^ '  which  has  been  used  in  assorted  form! 
to  study  attenuation  phenomena  with  a  fair  degree  of  succes; 
(Teng,  1968;  Solomon^^  1970;  and  several  others]. 

The  problem  is  simplified  to  a  point  source  in  an  eerU 
tor  which  geometric  ray  theory  and  linear  elasticity  and 
anelasti city  are  valid.  Then  the  observed  amplitude  spec¬ 
trum  of  an  isolated  body  wave  from  an  earthquake  may  be 
wr i tten : 


A(f)  =  S(f)  R  ( 9  ,  <J> )  A  (f)  A  (f) 


(1) 


where  S  is  the  amplitude  spectrum  at  the  source,  R  is 
f!..  radiation  pattern,  a  function  of  propagation  direction 
(.!),  Ap  is  the  transfer  function  for  propagation  through 
<urth  and  Ap  is  the  transfer  function  for  the  instru¬ 
ment.  Dy  geometric  ray  theory  „e  can  factor  A  (f)  as 
t  ol lows : 
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Ap(f)  =  G  •  exp[-ft*J  •  exp[-f6t*]  (,  ) 

where  G  is  a  frequency-independent  factor  which  accounts 
for  geometrical  spreading,  and  the  remainder  of  the  expres¬ 
sion  on  the  right-hand  side  of  (2)  accounts  for  attenuation. 

The  quantities  t*  and  6t*  are  given  by  integrals 
of  the  reciprocal  qualtiy  factor  over  the  ray  path.  We 
consider  that  Q  ^  at  any  point  within  the  earth  can  be 
written  as  a  sum  of  some  spherically  symmetric  model  Q~^(r) 
and  a  residual  term 

Q"1^)  =  Q~ 1  ( r )  +  5Q-1(r)  ,  (3) 

a  relationship  valid  for  small  losses.  In  general,  all 
quantities  in  (3)  may  also  depend  on  frequency.  We  then 
define 


t* 


v  1  (s) 


ds 


(4) 


and 


^L*  - 


-1 


-1 


f>Q  (s)  v  (s)  ds 


S 


(5) 
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where  ds  is  an  increment  of  length  along  the  ray  path 
S.  Equation  (2)  then  follows  from  the  definitions  (4)  and 
(5)  as  long  as  variations  in  v  ^  are  small  compared  to 


If  the  amplitude  spectrum  of  a  given  body  wave  is 
measured  at  two  stations,  then  the  ratio  of  the  spectra 
can  be  used  to  determine  the  difference  in  attenuation  of 
the  waves  to  the  two  respective  stations.  If 


R12(f) 


=  A1(f) 
A2(f) 


then 


In  R12(fj  +  f(t*  -  t*)  -  ln[I1(f)/I2(f)l  = 

c12  +  f(6t2  "  6tV  <6> 

The  Wt-lnnd  side  of  (6)  can  be  measured  or  calculated, 
and  should  vary  roughly  linearly  with  frequency  if  fit* 
is  a  slowl y-varyinq  function  of  frequency.  Thus  the  quanti¬ 
ty  (At*  -  At*  )  may  be  determined  by  fitting  a  straight 
to  the  lelt-hand  side  of  (6).  An  example  of  such  a 
measurement  i s  shown  in  Figure  3.  As  a  further  check, 
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thc  intercept  of  such  a  straight  line,  given  by  the  constant 
c12,  can  be  calculated  from  the  source  mechanism  of  the 
earthquake  and  a  suitable  velocity  model  for  the  earth. 

By  equation  (5) ,  the  contributions  to  6t*  can  arise 
from  any  portion  of  the  ray  path.  It  is  expected,  however, 
that  iSt*  will  be  controlled  primarily  by  the  upper  mantle, 
in  particular  the  as thenosphere ,  where  values  of  Q  1  are 
known  to  be  large  and  variable  [ Solomon ,  1972a]  and  where 
seismic  velocity  and  other  physical  properties  also  show 
lateral  variation.  This  intuitive  deduction  receives  support 
from  the  observation  that  station  differences  in  6t*  for 
shear  waves  from  deep-focus  earthquakes  do  not  generally 
depend  on  tie  source  location  [Solomon  and  Toksflz,  1970; 
Solomon ,  1971].  For  shallow  sources,  of  course,  the  body 
waves  to  stations  at  telescismic  distances  must  pass  through 
the  asthenosphere  twice.  Only  if  the  variation  in  upper 
mantle  attenuation  beneath  the  stations  is  known  a  priori 
or  if  the  source  is  located  in  a  region  of  unusually  low  Q 
can  the  near-source  and  near-receiver  contributions  to  ft* 
be  separated  with  some  confidence. 
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I) i  sens:,  i on  o T  Assumption:: 

.'IrVfi  .« I  <>l  I  In-  .  i :  t :  s  i  ini|  >  I  iniiM  iii.hIi-  iii  I  In-  ilcriv.lt  ion.-;  i>| 
tlu*  previous  section  either  are  incorrect  or  restrict  the 
permissible  data  set.  The  frequency  range,  .015  to  .17  Hz, 
of  the  shear  waves  considered  in  this  study  includes  fre 
quencies  too  low  for  geometric  ray  theory  strictly  to  apply 
and  too  high  for  the  earthquakes  indicated  above  to  be  re¬ 
garded  as  point  sources.  Further,  the  possibility  of  con¬ 
tamination  of  the  S  wave  by  other  phases  should  be  examined. 
We  discuss  below  the  degree  to  which  failure  of  these  assump¬ 
tions  may  affect  our  conclusions  about  seismic  attenuation. 

For  propagation  through  the  mantle,  use  of  geometric 
ray  theory  is  questionable  probably  only  if  the  wave  bottoms 
in  the  immediate  vicinity  of  a  transition  zone,  i.e.  a  sharp 
incrr  ase  in  velocity  with  depth.  In  this  work  we  consider 
S  waves  which  have  traveled  epicentral  distances  of  39  to 
69  degrees  or,  alternatively,  have  bottomed  in  the  mantle 
at  depths  of  about  900  to  1800  km.  Transition  zones  for 
shear-wave  velocity  at  such  depths  in  the  mantle  are  not 
particularly  we 11 -documented .  though  second  order  discontin- 
utitios  in  dt/dA  for  S  waves  have  been  suggested  to  exist  at 
A  ~  4J°  [Hales  and  Roberts,  1970a]  and  near  A  =  70°  [Fairborn, 
1969).  Presumably,  determinations  of  6t*  obtained  from 
waves  at  such  distances  should  be  used  with  care. 


The  amplitude  spectrum  of  a  body  wave  is  distorted  by 


passage  through  a  layered  crust.  This  distortion  can  be 
quite  large  for  periods  less  than  about  10  seconds;  it  is 
sensitive  to  crustal  structure;  and  it  tends  to  be  worse 
for  SV  waves  than  for  SH  waves  [Haskell ,  I960,  1962;  Ben- 
Menahem  et  al . ,  1965].  We  have  chosen  not  to  correct  shear 
wave  spectra  for  station  crustal  response,  primarily  because 
of  the  large  uncertainty  in  the  velocity  structure,  particu¬ 
larly  for  shear  waves,  in  the  crust  beneath  most  seismograph 
stations.  A  few  attempts  at  correcting  for  the  crust  con¬ 
vinced  us  that  the  uncertainty  in  the  correction  term  is 
comparable  to  the  correction  itself.  Further,  for  a  suffi¬ 
ciently  wide  frequency  band  the  effect  on  spectral  ratios 
of  different  station  crustal  responses  is  to  introduce 
several  peaks  and  troughs  without  appreciably  altering  the 
longer  term  trend  as  a  function  of  frequency. 

For  shallow  earthquakes,  the  body-wave  amplitude  spectrum 
will  be  multiplied  by  an  additional  transfer  function  because 
o t  the  crust  and  free  surface  near  the  source.  Earthquakes 
on  the  mid-ocean  ridge  system  are  known  to  be  very  shallow 
[Tsai,  1969;  Weidner ,  1972];  the  depth  of  the  2  June  1965 
event  was  constrained  by  Weidner  [1972]  to  be  only  3+2  km. 


Further, the  S  waves  reflected  from  the  free  surface  will  add 
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cons  true  tive  ly  to  the  direct  S  for  both  SV  and  SH  components 
for  the  two  earthquakes  considered  here.  Thus  the  ratio  of 
transfer  functions  associated  with  crustal  structure  at  the 
source  for  stations  at  two  somewhat  different  distances  will 
be  nearly  unity. 

A  mor  '  serious  complication  to  the  measurement  of  6t* 
is  the  finiteness  of  the  earthquake  source.  The  effects  on 
the  body-wave  radiation  pattern  of  rupture  propagation  on  a 
fault  plane  of  finite  extent  have  been  calculated  for  both 
the  frequency  and  time  domains  [Ben-Menahem,  1962;  Bollinger, 
1968;  Fukao,  1970].  These  effects  will  significantly  alter 
the  body-wave  amplitude  spectra  in  the  frequency  range  of 
interest  here  even  for  relatively  small  sources. 

As  a  guide  to  the  nature  of  the  error  in  <5t*  due  to 
neglect  of  finite  source  dimensions,  we  computed  the  varia¬ 
tion  with  azimuth  of  amplitude  spectrum  due  to  horizontal 
rupture  propagation  on  a  vertical  fault  [Ben-Menahen,  1962] 
and  the  resulting  dependence  on  azimuth  of  the  apparent 
differential  attenuation  obtained  from  the  slope  of  the 
ratio  of  the  amplitude  spectrum  at  the  given  azimuth  to 
the  spectrum  in  the  rupture  direction.  A  uniform  rupture 
velocity  of  2  km/sec  and  fault  lengths  of  10,  20  and  30  km 
were  used  a s  representative  values  for  mid-Atlantic  ridge 
earthquakes  of  magnitude  5.5  to  6  [e.g.  Udias ,  1971].  A 


constant  horizontal  phase  velocity  of  8  km/sec  was  used. 
This  corresponds  to  an  epicentral  distance  of  about  50°for 
S  waves  from  shallow  events;  the  range  in  phase  velocity 
for  the  distance  range  considered  in  this  study  is  about 
7.3  to  9.5  km/sec.  Amplitude  spectra  were  first  smoothed 
by  averaging  over  a  moving  window  .02  Hz  in  width  (see 
next  section) ,  and  apparent  differential  attenuation  was 
calculated  for  the  frequency  band  .01  to  .35  Hz.  The  value 
for  6t*  thus  determined  are  shown  in  Figure  4. 

The  effect  of  finiteness  is  to  introduce  troughs  in 
the  amplitude  spectra  at  frequencies  which  depend  on  azi¬ 
muth  and  on  the  difference  in  time  of  rupture  at  the  two 
ends  of  the  fault.  As  a  result,  there  will  be  paired 
peaks  and  troughs  in  the  spectral  ratio,  with  the  number 
of  such  features  dependent  on  the  fault  length  L  and  on 
the  frequency  band.  For  L  =  10  km,  the  lowest-frequency 
trough  is  located  at  f  =  .16  Hz  (4>  =  180°)  to  f  =  .27Hz 
=  °C)-  The  spectral  ratio  used  to  calculate  6t*  thus 
contains  a  single  trough  centered  at  a  frequency  above  the 
highest  frequency  in  the  band.  There  is  a  large,  azimuth- 
dependent  variation  in  the  apparent  differential  attenua¬ 
tion  as  a  result.  For  L  >  20  km,  there  is  at  least  one 
trough-peak  pair  in  each  spectral  ratio,  so  the  azimuthal 
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variation  of  St*  is  less  pronounced.  In  general,  the 
effect  of  source  finiteness  on  St*  decreases  as  fault 
length  increases  and  as  higher  frequencies  are  included 
in  the  frequency  band  [cf.  Ward  and  Toksflz,  1971].  The 
lowest  frequency  trough  is  in  the  observed  spectra  from 
the  event  of  13  February  1967  generally  occurs  at  a  fre¬ 
quency  in  the  range  of  .04  to  .10  Hz  and  is  azimuth- 
dependent,  though  there  is  considerable  scatter.  This  is 
consistent  with  a  fault  length  of  20  to  30  km,  and  indi¬ 
cates  that  the  effect  of  finite  source  size  on  St*  will 
probably  be  roughly  as  shown  in  Figure  4  for  these  values 
of  L.  For  the  2  June  1965  event,  the  geometry  is  some¬ 
what  different,  but  the  lowest-f requency  troughs  are  at 
about  the  same  or  slightly  higher  frequencies  than  for 
the  other  earthquake. 

It  thus  appears  likely  that  finite  source  dimensions 
and  rupture  velocity  will  contribute  a  reasonably  small 
(_  6  sec)  term,  which  varies  smoothly  with  azimuth,  to  the 
apparent  differential  attenuation.  If  observed  differences 
in  St*  are  much  larger  than  6  sec  or  if  there  are  large 
changes  in  St*  with  a  small  change  in  azimuth,  then  such 
observations  cannot  be  ascribed  to  source  effects  but  must 
reflect  genuine  variations  in  Q 
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The  S  waves  included  in  this  study  were  sufficiently 
well  isolated  from  other  body-wave  arrivals  (SS,  ScS ,  SKS  , 
etc.)  so  that  the  time  window  chosen  for  spectral  analysis 
does  not  include  these  phases.  A  further  potential  compli¬ 
cation  is  the  contamination  of  long-period  SV  waves  by  S- 
coupled  PL  waves  [Oliver ,  1961].  From  published  phase  and 
group  velocity  curves  for  fundamental-mode  PL  waves  [Chander 
et  al. ,  1968],  we  tested  whether  PL-wave  contamination  was 
likely  for  the  epicentral  distances,  frequency  bands,  and 
time  windows  employed  in  our  analysis  of  SV  waves.  With  a 
few  exceptions  noted  below,  all  involving  epicentral  dis¬ 
tances  less  than  50°  and  predominantly  oceanic  paths,  PL 
waves  should  not  have  affected  SV-wave  spectral  amplitudes. 


ATTENUATION  MEASUREMENTS 

Measurements  of  the  differential  attenuation  of  shear 
waves  were  made  in  the  same  fashion  as  in  Solomon  and  Toksflz 
11970].  Seismograms  from  horizontal  components  of  the  WWSSN 
wore  digitized  at  an  interval  of  about  0.7  sec,  rotated,  and 
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bandpass  filtered  after  mean  and  linear  trend  were  removed. 

SV  and  SH  components  wore  Fourier  transformed,  and  the  am¬ 
plitude  spectra  were  smoothed  by  applying  to  each  a  moving- 
average  window  of  width  .02  Hz.  The  smoothed  spectra  were 
corrected  for  the  instrumental  response,  using  the  theore¬ 
tical  formulas  of  Hagiwara  [1958].  Spectral  ratios  were 
calculated,  with  a  single  reference  station  for  each  event, 
and  were  corrected  for  geometric  spreading  using  travel¬ 
time  distance,  relationships  of  Hales  and  Roberts  [1970a] 
and  for  differences  in  fit*  (see  equation  6)  using  the  tech¬ 
nique  of  Julian  and  Anderson  [1968],  with  the  velocity  model 
SLUTD1  of  Hales  and  Roberts  [1970b]  and  a  Qq1  (r)  model  equal 
to  Tsai  and  Aki's  [1969]  modification  of  model  MM8  [Anderson 
et  al . ,  1965].  The  values  of  fit*  -  fit*  were  then  taken 
to  be  the  negative  of  the  slopes  of  the  straight  lines  fit 
by  least  squares  to  the  corrected  spectral  ratios. 

rn  Table  2  arc  given  the  measured  values  of  fit*  for 
shear  waves  from  the  t rans form- fault  earthquake  of  13  February 
1967.  Only  the  SH  component  was  used,  because  the  propaga¬ 
tion  direct  ions  of  interest  (within  a  few  tens  of  degrees 
azimuth  of  the  strike  of  the  fracture  zone)  are  near  a 
maximum  in  the  SH  radiation  pattern  but  near  a  nodal  plane 
for  SV.  The  list  of  stations  included  in  the  table  is 


exhaustive,  in  that  it  represents  all  WWSSN  stations  in  the 
distance  range  39 1  to  about  70°  that  recorded  clear  S  waves 
on  both  horizontal  components.  The  first  four  stations  in 
the  table  are  in  Europe  or  the  Middle  East,  the  next  six 
are  in  South  and  Central  America  oi  in  the  Caribbean,  and 
the  last  thirteen  are  in  North  America.  The  differential 
attenuation  fit*  .  is  determined  from  the  slope  of  the 
ratio  of  the  amplitude  spectrum  at  the  station  indicated 
to  that  at  the  reference  station  over  the  frequency  band 
given  at  the  bottom  of  the  table.  The  frequency  band  is 
as  wide  as  can  be  permitted  by  the  requirement  of  a  reason¬ 
able  signal-to-noise  ratio.  At  several  stations  with  the 
highest  values  of  6t*,  all  shear-wave  energy  at  frequencies 
above  .15  Hz  appears  to  have  been  removed  (e.g.  ALQ  in  Figure 
3) .  The  choice  of  reference  station  does  not  affect  the 
relative  values  of  fit*,  though  it  does  influence  the 
standard  deviation  of  the  straight-line  fits  to  the  spectral 
ratios.  These  standard  deviations,  also  given  in  the  table, 
are  only  lower  bounds  on  the  actual  uncertainty  in  (St*. 

Four  of  the  stations  lie  very  near  to  the  SH  nodal  plane  as 
predicted  from  Figure  2;  values  of  fit*  for  these  stations 
arc'  considerably  less  reliable  than  the  others. 
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Note,  from  equation  (6),  that  wo  can  determine  only 
differences  in  St*  between  stations.  A  constant  could 
be  added  to  all  entries  in  the  last  column  of  Table  2  with¬ 
out  altering  any  of  the  spectral  ratios.  The  stations  with 
larger  vaLues  of  <5t*,  corresponding  to  greater  attenuation 
would,  of  course,  remain  as  shown. 

As  a  rough  check  on  the  validity  of  our  interpretation 
of  spectra]  slopes,  we  compared  the  intercept,  c^2  in 
equation  ( f> )  ,  obtained  by  extrapolating  the  straight-line 
fit  to  i  -  0,  with  the  value  expected  from  the  radiation 
pattern  and  the  differences  in  geometrical  spreading  for 
the  paths  from  the  source  to  the  respective  pair  of  stations. 
Except  for  the  stations  near  a  nodal  plane,  which  all  re¬ 
corded  amplitudes  significantly  larger  than  predicted  from 
the  radiation  pattern,  all  values  of  c 12  were  within  a 
factor  of  ?  of  the  predicted  values  and  for  half  the  stations 
were  within  a  factor  of  1.5.  Considering  the  uncertainties 
in  the  corrections  and  the  omission  of  ground  motion  ampli¬ 
fication  by  local  geology  or  topography,  the  agreement  is 
fairly  good  and  lends  support  to  the  analysis  techniques 
adopted  in  this  work. 

Measured  vaLues  of  differential  attenuation  of  shear 
waves  for  the  ridge-crest  earthquake  of  2  June  1965  are  given 
m  Table  5.  For  this  event  only  the  SV  component  was  used; 
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SV  amplitudes  at  all  but  one  of  the  stations  J.n  Table  3 
are  greater  than  SH  amplitudes,  and  for  more  than  half  the 
stations  the  SV/SH  amplitude  ratio  is  3  or  more,  a  conse¬ 
quence  of  the  S  wave  radiation  pattern  from  a  normal-faulting 
earthquake.  As  with  the  previous  earthquake,  the  stations 
in  Table  3  include  all  usable  WWSSN  records  at  stations  in 
the  distance  range  40°  to  70°.  KTG  and  AKU  are  in  Green¬ 
land  and  Iceland,  respectively;  the  next  eight  are  in  Europe; 
SDB  1S  m  Africa;  the  next  three  are  in  South  America;  and 
the  last  nine  are  in  North  America.  The  time  window  and 
frequency  band  used  are  slightly  different  than  for  the 
previous  event,  because  of  both  relatively  more  energy  at 
higher  frequencies  and  the  need  to  minimize  interference 
from  the  PL  wave.  At  three  stations,  noted  in  the  table, 
there  may  still  be  contamination  from  PL  waves  at  the  lowest 
frequencies  in  the  passband. 

We  did  not  compare  the  intercepts  of  the  straight-line 
fits  to  the  spectral  ratios  with  their  expected  values  be¬ 
cause  no  fault-plane  solution  that  matched  P-wave  first 
motion  [e.g.  Sykes,  1970a]  could  fit  the  S-wave  polarization, 
specifically  the  direction  of  first  motion  of  the  SH  compo¬ 
nent..  The  problem  lies  in  the  standard  projections  used  to 
nuq>  surface  observations  of  first  motion  back  onto  the  lower 
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focal  hemisphere.  Because  of  pronounced  lateral  velocity 
gradients  beneath  the  ridge  crest,  body  waves  do  not  main¬ 
tain  a  constant  ray  parameter  along  their  path.  Such  lateral 
heterogeneity  can  also  account  for  the  apparent  non-orthogo¬ 
nality  of  nodal  planes  frequently  observed  for  normal  faultinq 
earthquakes  on  the  ridge  crest  (c.g.  Sykes ,  1967,  1970a]  and 
the  discrepancy  between  focal  mechanisms  obtained  from  P-wave 
first  motions  and  those  determined  from  surface  waves  [Weidner , 
1972].  This  subject  will  bo  treated  in  more  detail  elsewhere 
[Solomon  and  Julian,  1973]. 


DISCUSSION  AND  INTERPRETATION 

The  important  parameter  controlling  the  differential 
attenuation  measurements  presented  above,  at  least  controlling 
the  contribution  to  fit*  from  heterogeneous  Q  1  in  the 
upper  mantle  near  the  mid- Atlantic  ridge,  is  the  azimuth  at 
thf  source  of  the  path  to  each  station.  This  is  illustrated 
in  Figure  5,  whore  we  plot  fit*  versus  azimuth  for  the 
Gibbs  fracture  zone  earthquake  of  13  February  1967.  A 
constant  (12.4  sec)  has  been  subtracted  from  each  of  the 
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values  of  <St*  given  in  Table  2  so  that  the  mean  of  all 
values  is  zero. 

The  important  feature  of  Figure  5  is  the  band  of  very 
high  values  of  6t*  between  azimuths  of  273°  and  292°  ±  2° 
azimuth.  The  two  edges  of  this  band  are  interpreted  some¬ 
what  differently.  The  sharp  increase  in  <$t*  with  increasing 
azimuth  at  273°  azimuth  can  be  readily  explained  by  noting 
that  for  larger  azimuths  the  shear  waves  must  propagate  be¬ 
neath  the  crest  of  the  mid-Atlantic  ridge  while  for  lesser 
azimuths  the  waves  pass  under  only  older  sea  floor  (see 
Figure  1).  That  the  value  of  273°  azimuth  does  not  correspond 
precisely  with  the  fracture-zone  strike  or  the  local  direc¬ 
tion  of  plate  divergence  [Chase,  1972]  should  not  be  too 
disturbing  in  view  of  the  uncertainties  in  all  of  these  quan¬ 
tities.  The  conclusion  is  clear  that  Q_1  in  the  shallow 
mantle  beneat  i  the  mid-Atlantic  ridge  is  substantially  higher 
than  at  comparable  depths  beneath  sea-floor  25  to  30  million 
years  old. 

The  decrease  in  (St:*  with  increasing  azimuth  near 
2')2"  a/,  i  ninth  was  unexpected.  The  transition  is  sharp,  and 
t  In-  values  of  differential  attenuation  for  stations  at  azi¬ 
muths  of  204  to  330°,  even  though  S  waves  to  these  stations 
must  propagate  beneath  the  ridge  crest,  are  comparable  to 
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values  obtained  at  stations  in  the  northeast  and  southwest 
quadrants,  where  S  wave  paths  pass  beneath  only  relatively 
old  sea  floor.  We  ii fer  from  this  that  there  exists  beneath 
the  crest  of  the  mid-Atlantic  ridge  a  localized  zone  of 
very  low  Q.  The  boundaries  of  this  zone  must  be  very  sharp. 
The  shear  waves  leaving  the  fracture  zone  at  azimuths  in  the 
range  273°  to  292°  pass  through  this  low-Q  zone,  but  S  waves 
leaving  the  source  at  greater  azimuths  do  not. 

Some  reasonably  firm  limits  can  be  placed  on  the  dimen¬ 
sions  of  the  low-Q  zone.  As  a  working  hypothesis  we  will 
treat  this  zone  as  a  basically  two-dimensional  feature; 
i.e.,  one  that  is  more  or  less  continuous  beneath  the  ridge 
axis,  and  that  is  disrupted  only  by  offsets  in  the  ridge 
crest.  The  distance  between  the  ridge  crest  and  the  epi¬ 
center  is  an  upper  bound  on  the  half-width  of  the  low-Q 
zone,  since  clearly  the  earthquake  lies  outside  the  zone. 

The  ISC  and  CCS  epicenters  place  the  earthquake,  respectively, 
60  and  70  km  oast  of  the  ridge  crest  to  the  northwest.  Prob¬ 
ably  the  half-width  of  the  low-Q  zone  is  no  more  than  about 
60  km.  Further  the  low-Q  zone  must  be  shallow.  The  S-waves 
leaving  the  earthquake  at  azimuths  greater  than  292°  pass 
beneath  the  ridge  axis  at  depths  equal  to  or  greater  than  100 
to  ISO  km,  witli  the  range  in  these  depths  being  due  to  the 
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uncertainties  in  epicentral  coordinates,  source  depth,  and 
S-wave  velocity  model  near  the  ridge.  Thus  the  low-Q  zone 
immediately  beneath  the  ridge  crest  can  extend  no  deeper 
than  100  to  150  km  and  must  be  shallower  than  this  away 
from  the  crest. 

From  the  approximate  dimensions  of  the  low-Q  zone,  we 
can  estimate  a  lower  bound  on  0  ^  necessary  to  give  the  ob¬ 
served  difference  in  6t*  of  roughly  10  sec  between  S 

l)  II 

waves  that  did  and  did  not  pass  through  the  zone .  From 
equation  (5)  and  an  approximate  path  length  of  100  km  within 
the  low-Q  zone,  we  get  SQ  1  =  .1,  or  Q  =  10  within  the 
zone.  This  value  of  SQ  1  should  be  thought  of  as  only  a 
representative  value  for  the  frequency  band  used  in  this 
study,  in  view  of  the  likely  frequency  dependence  of  Q  in 
partially  melted  rock  [Walsh ,  1969;  Solomon ,  1972a].  Further, 
because  the  dimensions  of  the  low-Q  zone  are  only  upper 
bounds  and  because  the  highly  attenuated  shear  waves  may 
not  have  passed  through  the  full  length  of  the  zone,  Q  for 
S  waves  may  be  less  than  10  within  the  zone.  This  is  a 
very  low  value  for  Q,  sufficiently  low  that  the  usual  treat¬ 
ment  of  attenuation  as  a  small  perturbation  to  linear  elas¬ 
ticity  may  break  c  own . 

It  should  be  emphasized  that  the  high  values  for  6 t* 
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in  Figure  5  are  clearly  due  to  attenuation  in  the  upper 
mantle  near  the  source  and  cannot  be  explained  by  other 
source  or  path  effects.  From  Figure  4  and  the  discussion 
above,  it  may  be  readily  seen  that  rupture  propagation  on 
a  finite  fault  cannot  produce  an  apparent  differential 
attenuation  that  varies  as  rapidly  with  azimuth  as  does 
^t*  in  Figure  5.  Further,  lateral  variation  of  Q  1  in 
the  upper  mantle  beneath  the  receivers  cannot  explain  the 
observed  pattern  of  (St*.  The  stations  of  interest  for 
this  earthquake  are  all  in  North  Amer'ea,  where  the  attenu¬ 
ation  in  the  upper  mantle  has  been  well  studied  [Solomon 
and  Toksflz,  1970].  Correcting  the  values  of  6t*  in  Table 
2  for  the  relative  differences  in  upper-mantle  attenuation, 
where  known,  beneath  the  respective  stations  does  not  change 
by  a  substantial  amount  the  pattern  of  St*  versus  azimuth 
(see  Figure  ‘3)  .  Of  some  note  is  that  the  station  with  the 
highest  value  of  St*  after  those  corrections  is  FLO,  the 
nearest  station  to  the  epicenter  in  the  azimuth  band  270° 
to  290°  and  thus  the  station  observing  the  S  wave  with  the 
greatest  angle  of  incidence  at  the  source.  Other  bits  of 
evidence  supporting  the  low-Q  zone  are  the  shape  of  the  S- 
wave  spectra  for  waves  which  passed  through  the  zone  (the 
logarithm  of  the  spectral  amplitude  density  decreases  smoothly 
and  roughly  linearly  with  frequency,  as  in  Figure  3)  and  the 
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relatively  laige  fraction  of  emergent  arrivals  in  the  north¬ 
west  quadrant  of  the  fault-plane  solution  in  Figure  2. 

Not  all  of  the  shear-wave  spectra,  we  might  note,  de¬ 
creased  regularly  with  frequency  as  in  Figure  3.  At  stations 
in  the  northeast  and  southwest  quadrants,  the  spectra  appeared 
to  have  the  characteristic  "corner  frequency"  at  the  value 
(.07  Hz)  roughly  appropriate  to  source  dimensions  of  about 
30  km  [Brune ,  1970]  ,  consistent  with  the  considerations  of 
"directivity"  [Ben -Men ahem,  1962]  discussed  above. 

It  is  possible  that  properties  peculiar  to  the  Gibbs 
fracture  zone,  particularly  if  it  is  a  leaky  transform  fault, 
may  contribute  to  the  pattern  of  attenuation  given  in  Figure 
5.  The  high  value  of  6t*  at  ATH  would  support  such  a 
notion,  though  this  may  just  as  reasonably  be  a  near— receiver 
effect.  The  extent  to  which  a  shallow  zone  of  very  low  Q 
is  a  universal  feature  of  mid-oceanic  ridges  can  be  ascer¬ 
tained  only  by  further  observations. 

The  measured  values  of  6t*sv  for  the  ridge-crest 
earthquake  of  2  June  1965  are  plotted  as  a  function  of  azimuth 
in  Figure  6.  Again,  a  constant  value  (0.9  sec)  has  been 
added  before  plotting  to  all  of  the  values  of  6t*  in  Table 
\  so  that  the  average  value  of  fit*  for  this  data  set  is 
/.o  ro . 

In  at  least  the  two  northern  quadrants  of  Figure  6, 
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there  is  a  small  azimuthal  dependence  for  6t*.  The  dif¬ 
ferential  attenuation  is  higher  by  about  5  sec,  on  the 
average,  for  propagation  paths  more  nearly  parallel  to  the 
ridge  axis  (i.e.  to  European  stations)  than  for  paths  per¬ 
pendicular  to  the  strike  of  the  ridge  (to  North  American 
stations).  This  could  be  explained  if  Q  in  the  asthenosphere 
increases,  or  if  a  low-Q  zone  decreases  in  thickness,  as  a 
function  of  distance  from  the  ridge  crest. 

For  this  earthquake,  we  cannot  rule  out  explanations 
other  than  true  anelastic  losses  for  the  pattern  of  ^t* 
in  Figure  6.  From  Figure  4,  for  instance,  the  possibility 
that  fit*  is  controlled  by  the  soiree  must  be  admitted. 
Further,  the  difference  in  the  distributions  of  Q 
between  the  upper  mantles  of  Europe  and  North  America  is 
not  known.  This  probably  will  not  provide  an  explanation 
for  Figure  6,  however,  because  most  of  the  North  American 
stations  used  are  in  western  United  States,  when  upper 
mantle  Q  is  anomalously  low  [Solomon  and  Toksflz,  1970; 
Solomon,  1972a].  because  Q  in  the  upper  mantle  beneath 
most  of  Europe  is  likely  at  least  as  large  as  that  under 
western  U.S.,  the  difference  in  fit*  between  stations  in 
the  northeast  and  northwest  quadrants  would  be  accentuated 
if  corrections  for  near-receiver  attenuation  were  applied 


to  6t*  values  at  the  Kuropean  stations.  A  final  compli¬ 
cation  is  that  the  few  data  in  the  southern  two  quadrants 
do  not  mirror  the  pattern  of  the  northern  two,  though  this 
can  be  reasonably  ascribed  to  the  offset  in  the  ridge  crest 
south  of  the  epicenter  as  noted  earlier. 

As  an  aside,  it  might  be  pointed  out  that  three  of  the 
determinations  of  <5 1 *  that  appear  most  anomalous  when  com¬ 
pared  with  values  at  adjacent  azimuths  (ATM  and  SJG  in  Table 
2  and  Figure  5,  MAL  in  Table  3  and  Figure  6)  are  at  stations 
about  42°  from  the  epicenter.  As  noted  earlier  this  epicen- 
tral  distance  may  mark  a  discontinuity  in  the  slope  of  the 
travel  time  curve,  so  that  geometric  ray  theory  may  be 
inappropriate . 

The  main  conclusion  to  be  derived  from  Figure  6  is  that 
there  are  no  large  variations  with  propagation  direction  of 
attenuation  of  shear  waves  from  an  earthquake  on  the  ridge 
crest.  Since  our  measurements  reflect  only  differences  in 
total  attenuation.  Figure  6  is  still  consistent  with  the 
evidence  discussed  above  for  a  zone  of  very  low  Q  beneath 
the  ridge  axis  as  long  as  all  of  the  shear  waves  considered 
tor  this  event  passed  through  the  zone.  This  requirement 
allows  us  to  place  a  lower  bound  on  the  thickness  of  the 
•/.one.  The  bound  will  depend  on  the  depth  and  on  the  velo¬ 
city  model  near  the  source.  For  the  stations  closest  to  the 
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op  i  center ,  the  angle  of  incidence  in  the  shallow  mantle  is 
30°  to  40°  for  any  reasonable  upper-mantle  velocity  model. 
This  implies  the  half-width  of  the  low-Q  zone  must  be  at 
least  50  to  60  km  at  100  km  depth  or  at  least  25  to  30  km 
at  half  that  depth.  Since  these  values  are  comparable  or 
only  somewhat  smaller  than  the  upper  bounds  on  the  half¬ 
width  of  the  low-Q  zone  derived  earlier,  the  boundaries  of 
the  zone  are  fairly  well  constrained. 

What  phenomenon  will  give  rise  to  a  sharply  defined 
zone  of  extremely  low  Q?  The  simplest  explanation  is  that 
the  high  attenuation  is  associated  with  partial  melting, 
as  discussed  in  the  introduction,  and  that  the  low-Q  zone 
corresponds  to  a  region  where  the  melt  phase  occupies  a 
relatively  large  fraction  of  the  volume.  A  partially  melted 
region  beneath  ridge  crests  has  previously  been  predicted 
from  models  of  convective  upwelling  beneath  the  ridge  [Bott , 
1965;  Qxburgh  and  Turcotte ,  1968;  Wyllie,  1971;  Forsyth  and 
P_ross,  1971].  Wyllie  [1967],  in  particular,  emphasized  the 
distinction  between  incipient  melting,  when  temperatures 
lie  above  the  solidus  of  mantle  material  partially  saturated 
with  water  but  below  the  anhydrous  solidus  and  onlv  a  rela¬ 
tively  small  melt  concentration  is  present,  and  "normal"  or 
"dry"  melting,  when  temperatures  exceed  the  anhydrous  solidus 
of  mantle  material  and  large  quantities  of  melt  are  present. 
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Froin  Wyllie 's  discussion  and  the  extremely  low  Q  present 
beneath  the  ridge  crest,  we  identify  the  low-Q  zone  with 
a  region  of  "dry"  melting  and  the  boundaries  of  the  zone 
with  the  anhydrous  solidus  of  mantle  material. 

Such  an  explanation  is  very  plausible  from  the  stand¬ 
point  of  current  models  for  the  temperature  field  in  the 
mantle  beneath  mid-ocean  ridges.  In  Figure  7  is  shown 
such  a  temperature  field,  taken  from  a  convection  model 
of  Andrews  [1972]  incorporating  variable  viscosity  and  a 
lithosphere-spreading  half-rate  of  1.2  cm/yr.  From  this 
temperature  field  and  from  the  schematic  phase  diagram  of 
Wyllie  [1971]  for  peridotite  in  the  presence  of  .1  percent 
water,  the  predicted  regions  of  incipient  and  "dry"  melting 
are  as  shown  in  Figure  7.  Similar  curves  were  also  given 
by  Wyllie  and  by  Oxburgh  and  Turcotte  [1968],  both  using 
the  temperature  field  of  the  latter  authors.  The  broad 
region  of  incipient  melting  corresponds  to  the  normal  oce¬ 
anic  asthonosphere ,  and  the  region  of  "dry"  or  "normal" 
melting  is  confined  to  depths  less  than  100  km. 

Since  the  model  of  Andrews  [1972]  incorporated  some 
arbitrary  assumptions  and  simplifications  and  Wyllie's 
| 1971]  i sopleth  was  not  intended  to  be  precise,  it  is  not 
very  meaningful  to  derive  quantitative  conclusions  from 
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Figure  7.  Nonetheless  it  is  apparent  that  the  region  of 
large  melt  concentration  is  much  wider  in  Figure  7  than  the 
low-Q  zone  inferred  above.  This  is  also  a  feature  of  the 
diagrams  of  Oxburgh  and  Turcotte  [1968]  and  of  Wyllie  [1971], 
Andrews  [personal  communication,  1972]  has  noted  that  the 
temperature  field  scales  approximately  as  the  thermal  gra¬ 
dient  in  the  lithosphere,  which  is  governed  by  the  somewhat 
arbitrary  values  for  normal  oceanic  heat  flux  and  for  thermal 
conductivity  in  the  lithosphere.  For  instance,  if  the  con¬ 
ductivity  in  the  lithosphere  is  raised  by  6  percent  from  the 
value  used  by  Andrews,  then  the  zone  of  "dry"  melting  shrinks 
to  the  shaded  region  of  Figure  7,  50  km  in  half-width  and 
shallower  than  60  km  in  depth.  Though  the  match  of  the 
dimensions  of  this  region  to  those  estimated  for  the  zone 
of  very  low  Q  may  be  accidental,  it  nonetheless  illustrates 
how  meauremonts  of  seismic  attenuation  can  serve  as  useful 
constraints  on  models  of  mantle  flow  and  temperature. 
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OTHER  RELATED  OBSERVATIONS 

There  have  been  a  number  of  studies  and  observations 
which  relate  to  the  findings  of  this  paper.  All  of  them 
are  more  or  less  consistent  with  a  shallow,  narrow  zone  of 
low  Q  associated  with  partial  melting  at  temperatures  in 
excess  of  the  dry  solidus. 

Molnar  and  Oliver  [1969]  noted  that  Sn  waves  which  have 

propagated  across  mid-ocean  ridge  crests  are  always  highly 

attenuated;  but  for  S  waves  from  fracture-zone  earthquakes, 

n 

paths  not  crossing  the  ridge  often  show  efficient  trans¬ 
mission.  Their  findings  differed  somewhat  from  ours  in 
that  some  transform-fault  earthquakes  (including  two  on 
the  Gibbs  fracture  zone)  located  further  from  the  ridge 

crest  than  60  to  70  km  displayed  inefficient  s  trans- 

n 

mission  for  paths  without  ridge  crossings.  Su'h  trans¬ 
mission  studies  may  not  be  directly  correlatable  with  our 
work,  however,  for  at  least  two  reasons.  The  higher  fre¬ 
quencies  of  Sn  waves  means  that  lower  values  of  Q  1  are 
needed  to  damp  out  most  of  the  energy;  thus  Molnar  and 
Oliver's  observations  may  bear  more  on  the  general  eleva¬ 
tion  of  the  asthenosphere  beneath  mid-ocean  ridges  (see 
K inure  7)  than  on  a  more  localized  region  of  exceedingly 
low  n.  Kurt  her,  the  apparent  attenuation  of  Sn  waves  can 


-108- 


probably  be  markedly  increased  by  negative  velocity  gradients 
with  depth,  associated  with  high  thermal  gradients,  in  the 

uppermost  mantle  [Hill ,  1971]. 

From  the  apparent  velocities  of  P  waves  recorded  in 
Iceland  from  earthquakes  on  the  mid-Atlantic  ridge,  Francis 
[1969]  deduced  that  anomalously  low  velocities  extend  to 
250  km  depth  in  the  mantle  within  150  km  of  the  ridge  axis. 
Such  a  half-width  and  depth  are  greater  than  for  the  low-Q 
zone  mentioned  above,  though  this  P-wave  data  may  also  pri¬ 
marily  reflect  lateral  variations  within  the  deeper  astheno- 
sphere  or  may  merely  indicate  that  the  upper  mantle  beneath 
Iceland  is  different  from  other  portions  of  the  mid-Atlantic 
ridge . 

The  differential  attenuation  of  shear  waves  from  an 
earthquake  on  the  crest  of  the  mid-Arctic  (Nansen)  ridge  was 
reported  by  Solomon  and  Toksflz  [1970],  though  the  number  of 
data  and  the  azimuthal  coverage  were  more  limited  than  in 
the  present  study.  After  correcting  for  differences  in 
upper-mantle  attenuation  beneath  the  receivers,  the  values 
of  At*  show  a  marked  azimuthal  dependence  (Figure  8), 
possLbly  indicating  a  difference  in  Q  1  between  the  upper 
mantle  beneath  the  oceanic  segment  of  the  ridge  and  that 
beneath  the  continental  shelf  or  else  unusually  high  Q  at 
several  hundred  kilometers  depth  beneath  the  Lomonosov  ridge. 
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Ward  and  Toksflz  [1971]  made  a  pair  of  observations 
similar  to  those  reported  in  this  paper.  They  determined 
the  difference  in  attenuation  of  P  waves  from  two  mid-Atlantic 
earthquakes  recorded  at  the  two  seismic  arrays  LASA  and 
NORSAR.  For  an  event  west  of  the  ridge  crest,  the  P  wave 
which  propagated  beneath  the  ridge  (to  NORSAR)  was  more 
attenuated;  while  for  a  ridge-crest  earthquake,  the  difference 
in  attenuation  at  the  two  stations  was  not  statistically  sig¬ 
nificant.  The  former  result  does  not  bear  directly  on  the 
proposed  low-Q  zone  beneath  the  ridge,  because  the  P  wave 
to  NORSAR  passed  under  the  ridge  crest  at  about  450  km  depth. 

It  has  generally  been  noted  that  body-wave  magnitudes 
of  ridge  crest  earthquakes  are  often  anomalously  lower  than 
values  for  earthquakes  of  comparable  surface-wave  magnitude 
located  elsewhere  [e.g.  Solomon ,  1972b].  Frequently,  in  fact, 
surface  waves  are  observed  from  mid-ocean  ridge  earthquakes 
for  which  no  body  waves  are  detected  [Marshall ,  1970;  Sykes , 
1970b] .  This  phenomenon  can  readily  be  explained  by  a  low- 
Q  zone  in  the  shallow  mantle  beneath  ridges;  the  short- 
period  body  waves  which  pass  through  the  zone  are  sharply 
attenuated  but  surface  waves  generated  by  the  same  earth¬ 
quake  are  little  affected.  Francis  and  Porter  [1972] 
recent  ly  reported  some  related  observations.  Two  ocean- 
bottom  seismometers  30  km  apart,  one  in  the  median  valley 


of  the  mid-Atlantic  ridge  and  one  off  the  crest,  recorded 
markedly  different  rates  of  microearthquake  seismicity.  The 
much  lower  rate  of  activity  off  the  crest  can  be  attributed 
to  low  Q  immediately  beneath  the  ridge  axis.  Francis  and 
Porter  also  invoked  low  Q  to  explain  why  their  seismometers 
did  not  record  arrivals  from  several  large  earthquakes  which 
occurred  at  relatively  close  range  during  the  recording 
period . 

Kay  et  al.  [1970]  have  presented  chemical  evidence 
indicating  that  rocks  at  mid-ocean  ridges  are  the  product 
of  extensive  partial  melting  at  shallow  depths  beneath  ridge 
crests.  They  suggest  that  30  percent  melting  at  depths  of 
15  to  25  km  is  indicated.  Certainly  such  large  melt  concen¬ 
trations  are  consistent  with  and  perhaps  required  by  a  Q  for 
shear  waves  of  less  than  10,  though  most  theoretical  models 
of  attenuation  in  partially  melted  rock  [e.g.  Walsh,  1969] 
break  down  for  such  large  volume  concentrations  of  the  liquid 
phase.  The  molting  model  of  Kay  et  al.  has  extensive  melting 
confined  to  a  narrow  dike  (less  than  10  km  wide) ,  in  disa¬ 
greement  with  our  observation  that  differential  attenuation 
of  shear  waves  from  a  ridge-crest  earthquake  shows  no  large 
dependence  on  propagation  direction.  The  chemical  data, 
however,  do  not  preclude  as  much  as  10  percent  melting  in 
the  upper  mantle  at  considerable  distance  from  the  ridge 
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crest  [Kay  et  al . ,  1970];  such  melt  concentrations  would 
still  give  rise,  in  all  probability,  to  very  low  Q.  That 
extensive  melting  may  persist  no  deeper  than  25  km  beneath 
the  ridge  axis,  as  proposed  by  Kay  et  al .  ,  is  not  precluded 
by  our  attenuation  measurements,  though  for  the  bottom  of 
the  low-Q  zone  to  be  so  shallow  would  favor  a  larger  value 
for  the  width  of  the  zone,  about  100  km,  and  would  require 
even  lower  values  of  Q  than  discussed  above. 


CONCLUSIONS 


There  are  pronounced  lateral  inhomogeneities  in  seismic 
attenuation  beneath  mid-ocean  ridges.  The  variation  with 
propagation  direction  of  the  attenuation  of  shear  waves 
from  an  earthquake  on  the  Gibbs  fracture  zone  indicates  the 
existence  of  a  zone  of  very  low  Q,  a  zone  no  wider  than  about 
100  km  and  no  deeper  than  100  to  150  km,  beneath  the  axis  of 
the  mid-Atlantic  ridge.  That  such  a  conclusion  could  be 
made  is  due  to  a  fortuitous  combination  of  features  of  this 
particular  earthquake:  (i)  the  event  was  offset  from  the 
ridge  crest;  (ii)  the  azimuthal  directions  from  the  source 
•o  (lie  boundaries  of  the  low-Q  zone  corresponded  to  a  lobe 
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in  the  SH  radiation  pattern;  (iii)  this  same  range  of 
azimuths  included  S  wave  paths  to  a  relatively  dense  net¬ 
work  of  stations  at  appropriate  distances;  (iv)  the  upper- 
mantle  attenuation  beneath  the  stations  was  known.  Wo  might 
note  that  neither  P  nor  SV  waves  from  transform-fault  earth¬ 
quakes  can  be  used  to  delineate  the  low-Q  zone,  because  the 
direction  along  the  fracture  zone  corresponds  approximately 
to  both  P  and  SV  nodal  pianos.  There  have  been  at  most  a 
handful  of  earthquakes  on  the  mid-ocean  ridge  system  in 
the  last  decade  which  satisfy  all  of  the  features  listed. 

Such  earthquakes  merit  considerable  further  study. 

We  observed  little  variation  with  propagation  direction 
of  attenuation  of  shear  waves  from  an  earthquake  on  the  crest 
of  the  mid-Atlantic  ridge ,  indicating  that  all  waves  have 
seen  a  similar  Q  environment.  If  the  earthquake  occurred 
over  a  shallow  low-Q  zone,  a  hypothesis  strengthened  by 
other  anomalous  properties  of  body-wave  amplitudes  from 
ridge-crest  events,  then  such  a  zone  is  wider  than  about  50 
km.  Thus  the  combination  of  observations  from  the  two  earth¬ 
quakes  allows  us  to  place  reasonably  tight  constraints  on  the 
dimensions  of  the  low-Q  zone. 

The*  very  low  values  of  Q  (10  or  less)  within  such  a 
zone  can  most  readily  be  explained  as  due  to  extensive  par¬ 
tial  melting,  probably  because  convective  upwelling  beneath 
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the  ridge  axis  has  raised  the  peridotitic  material  of  the 
upper  mantle  above  its  anhydrous  solidus.  Thus  the  boundaries 
of  such  a  low-Q  zone,  when  combined  with  phase  diagrams  of 
likely  mantle  materials,  can  serve  as  quantitative  input 
to  refine  models  of  temperature  and  flow  beneath  mid-ocean 
ridges . 

The  Q  model  proposed  in  this  paper  for  the  shallow 
mantle  beneath  radges  can  readily  be  tested  by  further  ex¬ 
periments,  the  most  fruitful  of  which  are  likely  to  be  short- 
range  seismic  propagation  studies  near  spreading  centers 
using  either  natural  or  artificial  sources  and  ocean-bottom 


sensors . 
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Table  2. 

Differential 

attenuation  of  S 

waves  from  the 

earthquake  of 

13  February  1967 

• 

Station 

Azimuth 

Distance 

6t* 

SH 

degrees 

degrees 

sec 

MSH 

64 

64.2 

3.0  ±  1.63 

TAB 

71 

55.5 

6.5  ±  1.0 

JER 

84 

53.3 

9.1  ±  1.1 

ATH 

87 

42.2 

20.2  ±  1.8 

TRN 

218 

47.4 

12.4  ±  1.0 

CAR 

224b 

49.7 

3.6  ±  1.2 

SJG 

228b 

42.6 

-0.5  ±  1.4 

BOG 

229b 

58.1 

10.4  ±  1.2 

QUI  *. 

231b 

64.4 

10.5  ±  1.9 

BHP 

237 

57.0 

9.5  ±  1.3 

ATL 

262 

40.4 

12.3  ±  2.4 

OXF 

267 

42.7 

10.5  ±  2.3 

JCT 

272 

51.8 

14.1  ±  1.3 

FLO 

273 

40.5 

25.0  ±  1.5 

LUB 

276 

50.9 

19.2  ±  2.3 

ALQ 

281 

52.7 

22.1  +  1.2 

TUC 

281 

57.1 

21.4  ±  2.6 

DUG 

290 

52.9 

27.6  +  0.9 

BOZ 

294 

48.8 

12.4  ±  2.0 

COR 

300 

55.7 

9.1  ±  2.1 

LON 

301 

53.5 

7.9  ±  2.1 

CMCC 

324 

39.5 

0 

COL 

330 

52.0 

9.1  ±  1.8 

standnrd  deviation  of  straight-line  fit  to  spectral  ratio 
mar  Sll  nodal  plane  Time  window  length:  58  ±  1  sec 

Filter  bandpass:  -015  -  .15  Hz 


c 


reference  station 


|  ■» 
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Table  3.  Differential  attenuation  of  S  waves  from  the 
earthquake  of  2  June  1965. 


Station 

Azimuth 

degrees 

KTG 

10 

AKU 

14 

KONb 

29 

ESK 

30 

NUR 

31 

VALC 

31 

STU 

41 

TRI 

46 

TOLC 

49 

MALC 

53 

ATH 

55 

SDB 

114 

PEL 

205 

ANT 

211 

NNA 

229 

JCT 

296 

TUC 

298 

DAL 

300 

ALQ 

302 

OXF 

304 

BKS 

305 

BOZ 

314 

MDS 

316 

AAM 

317 

Distance 

6t* 

sv 

degrees 

sec 

56.5 

-1.1 

± 

1.4a 

53.3 

5.8 

± 

2.1 

59.6 

0 

51.5 

1.0 

± 

1.9 

67.1 

-2.2 

± 

1.5 

46.1 

4.0 

+ 

1.6 

55.9 

7.1 

+ 

1.6 

58.2 

5.1 

± 

2.0 

44.0 

3.8 

± 

1.8 

42.8 

-10.9 

± 

2.1 

65.0 

1.5 

± 

2.6 

67.0 

8.6 

± 

2.5 

53.9 

-7.7 

± 

1.7 

45.7 

-4.3 

± 

2.0 

40.7 

-3.8 

± 

2.2 

50.5 

-2.6 

± 

1.3 

59.9 

-1.0 

± 

1.6 

48.2 

-7.1 

± 

1.7 

56.4 

-2.3 

+ 

0.9 

42.5 

-6.9 

+ 

3.8 

69.1 

-0.7 

4; 

1.5 

61.3 

3.7 

± 

1.4 

45.6 

-3.4 

± 

2.6 

41.1 

-9.2 

+ 

2.5 

a  standard  deviation  of  straight-line  fit  to  spectral  ratio 
b  reference  station 

r  possible  interference  with  oceanic  PL  waves 


Time  window  length: 


54  *  1  sec;  Filter  bandpass:  .028 


175  Hz 
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FIGURE  CAPTIONS 

Fig.  1.  Location  of  earthquakes  considered  in  this  study 
(large  dots) .  Small  dots  are  other  epicenters  on 
the  mid-Atlantic  ridge  located  by  USCGS/NOAA  for 
the  years  1962-69.  The  small  arrows  denote  the 
sense  of  relative  movement  of  the  adjoining  plates 
for  the  event  of  13  February  1967  and  the  inferred 
axis  of  maximum  tension  for  the  event  of  2  June  1965 
[Sykes  ,  19  70a]  . 

Fig.  2.  Fault  plane  solution  from  P-wave  first  motions  for 
the  earthquake  of  13  February  1967  (equal-area  pro¬ 
jection)  .  Open  circles  are  dilatations,  closed 
circles  are  compressions,  and  crosses  indicate 
proximity  to  a  nodal  plane.  Smaller  symbols  are 
somewhat  less  reliable  than  the  larger  ones,  due  to 
lower  signal -to- noise  ratio  or  the  emergent  nature 
of  the  arrival.  All  data  were  read  from  long-period 
records  of  stations  in  the  WWSSN  or  the  Canadian  net¬ 
work.  The  strike  <J>  and  dip  6  of  the  nodal  planes 
are  also  given.  If  the  conversion  from  distance  to 
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angle  of  incidence  [Herrin,  1968;  a  velocity  of 
8  km/sec  was  assumed  at  the  source]  is  correct, 
then  4>  and  6  are  determined  to  within  about 
one  degree;  the  actual  uncertainty  is  probably 
larger . 

Fig.  3.  A  representative  determination  of  6t*^.  Shown 

are  the  SH-wave  amplitude  spectra  at  stations  ALQ 
and  CMC  (reference  station) ,  corrected  for  instru¬ 
ment  response.  Note  the  semi-log  scale.  For 
clarity,  the  ALQ  spectrum  (right-hand  scale)  has 
been  shifted  downward  by  a  factor  of  10  with  respect 
to  the  CMC  spectrum  (upper  left  scale) .  The  spectral 
ratio,  corrected  for  the  difference  in  epicentral 
distance  between  CMC  and  ALQ  (including  attenuation 
and  geometric  spreading)  and  for  the  source  radi¬ 
ation  pattern,  is  shown  over  a  more  limited  fre¬ 
quency  band.  The  differential  attenuation  is 
determined  by  a  least-squares  fit  of  a  straight 
line  to  the  spectral  ratio;  6t*  is  the  negative 
of  the  slope. 

Fig.  4.  The  effect  of  source  finiteness  on  measurements 

S’ 


of  <St* 


The  central  diagram  is  in  polar 
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coordinates  (r,4>)  ,  with  r  given  by  6t*  (see 
radial  scale,  in  sec)  and  <J>  equal  to  the  differ¬ 
ence  in  azimuth  between  the  rupture  direction  and 
the  direction  of  shear-wave  propagation.  Ampli¬ 
tude  spectra  were  calculated  following  Ben  Menahem's 
[1962]  relations  for  a  vertical  fault  with  horizontal 
rupture  propagation,  assuming  flat  source  spectra 
for  the  equivalent  point  source  and  using  a  rupture 
velocity  of  2  km/sec,  an  S-wave  horizontal  phase 
velocity  of  8  km/sec,  and  the  fault  length  L  as 
indicated.  Spectra,  smoothed  by  averaging  over  a 
.02  Hz  window,  are  shown  in  the  insert  for  the 
three  values  of  L  at  4>  =  0 .  The  differential 
attenuation  <5t*s(4>)  is  calculated  from  the  slope 
of  the  ratio  of  the  amplitude  spectrum  at  azimuth 
<t>  to  that  at  4>  =  0  over  the  frequency  band 
.01  to  .15  Hz. 

Fig.  5.  Differential  attenuation  6t*  of  SH  waves  from 
the  Gibbs  fracture  zone  earthquake.  The  plot  is 
in  polar  coordinates  (r , )  ,  with  r  given  by 
6t*  (see  radial  scale,  in  sec)  and  <j)  equal  to 
the  azimuth  of  the  respective  ray  path  at  the  source. 
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Closed  circles  are  the  measured  values  from  Table 
2,  with  a  uniform  constant  subtracted  to  give  the 
data  set  zero  mean.  The  error  bars  are  the  standard 
deviations  from  Table  2.  Open  circles  are  the 
measured  values  of  6t*  corrected  for  differences 
in  S-wave  attenuation  in  the  upper  mantle  beneath 
the  receivers ,  where  known  [Solomon  and  Toksflz, 

1970] .  The  strike  of  the  fracture  zone  [Fleming 
et  al . ,  1970],  identical  to  the  strike  of  the  fault 
plane  for  this  earthquake,  and  the  local  direction 
of  plate  divergence  [Chase ,  1972]  are  indicated  by 
arrows . 

Fig.  6.  Differential  attenuation  6t*  of  SV  waves  from 
the  earthquake  of  2  June  1965.  Coordinates  and 
symbols  are  as  in  Figure  5.  The  approximate  strike 
of  the  ridge  axis,  taken  from  earthquake  epicenters, 
and  the  local  direction  of  plate  divergence  [Chase , 
1972]  are  indicated  by  arrows. 

Fig.  7.  The  shape  of  a  zone  of  extensive  melting  beneath 
a  mid-ocean  ridge.  The  regions  of  incipient 
molting  (temperatures  between  the  water-saturated 
solidus  and  the  anhydrous  solidus)  and  of  "dry" 
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melting  (temperatures  between  the  anhydrous  solidus 
and  the  liquidus)  are  determined  from  the  temper¬ 
ature-depth  field  in  a  mantle-convection  model 
[Andrews ,  1972]  designed  to  simulate  sea-floor 
spreading  on  the  mid-Atlantic  ridge  and  the  schematic 
phase  diagram  of  Wyllie  [1971]  for  peridotite  in 
the  presence  of  .1  percent  water.  The  stippled 
area  is  the  region  of  "dry"  melting  if  the  thermal 
conductivity  in  the  lithosphere  is  6  percent  greater 
than  that  used  by  Andrews  (see  text)  . 

Fig.  8.  Lateral  variation  of  near-source  attenuation  for 

shear  waves  from  the  Arctic  earthquake  of  25  August 
1964.  The  difference  between  6t*g  determined  at 
several  North  American  stations  and  6t*  (S.  Am.) 

determined  from  deep  South  American  earthquakes 
recorded  as  the  same  stations  is  shown  as  a  function 
of  azimuth  at  the  source  at  right.  All  data  are 
from  Solomon  and  Toksfiz  [1970] .  The  projections 
onto  the  earth's  surface  of  several  S-wave  propa¬ 
gation  paths  to  North  America  are  indicated  on  a 
map  of  the  Arctic  region  (Lambert  equal-area 
projection)  at  left.  The  Lomonosov  ridge  is  well- 
defined  by  the  2-km  isobath  [simplified  from  Vogt 
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and  Ostenso,  1970]  while  the  deeper  mid-Arctic 
(Nansen)  ridge  follows  the  trend  of  earthquake 
epicenters  (crosses;  USCGS  and  NOAA  determina¬ 
tions)  . 
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3 . 7  Deep  Structure  and  Geophysical  Processes  Beneath 

Island  Arcs  by  N.H.  Sleep  (Abstract) 

The  deep  structure  and  geophysical  processes  associated 
with  island  arcs  were  examined  as  a  problem  in  heut  and 
mass  transfer.  Particular  attertion  was  given  to  the 
thermal  behavior  of  lithosphere  as  it  descends  into 
mantle,  seismic  transmission  through  descending  slabs, 
the  origin  of  the  magmas  which  erupt  on  island  arcs, 
and  the  cause  of  crustal  spreading  in  intra-arc  basins. 

The  factors  affecting  the  thermal  behavior  of  a 
lithosphere  slab  descending  into  the  mantle  are  so  com¬ 
plicated  and  numerous  that  only  numerical  methods  can 
accurately  account  for  them.  A  series  of  finite-differ¬ 
ence  calculations  of  the  temperature  field  was  made  to 
test  the  affect  of  phase  changes,  the  dip  of  the  slab, 
the  thermal  conductivity,  and  the  initial  geotherm. 
Reasonable  variations  in  those  parameters  and  geolo¬ 
gically  permissible  variations  in  radioactive  heating 
and  adibatic  compression  did  not  greatly  affect  the 
gross  thermal  structure  calculated  for  the  slab.  Ir¬ 
regularities  in  slab  movement  on  a  time  scale  of  less 
than  2  million  years  do  not  significantly  affect  the 
thermal  field.  Possible  factors  limiting  the  maximum 
depth  of  seismic  zones  include  thermal  assimilation  of 
the  slab  due  to  latent  heat  from  phase  changes  and  thermal 
conduction,  mechanical  disruption  of  the  slab,  and 
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the  lack  of  stress  at  great  depths  due  to  lower  density 
contrast  between  the  slab  and  the  mantle. 

Theoretical  ray  paths  through  the  numerically  cal¬ 
culated  thermal  models  of  slabs  were  computed.  The 
results  were  in  good  agreement  with  observed  travel  times. 
First  motion  amplitudes  of  P-waves  at  teleseismic 
distances  were  measured  from  long  and  short  period 
WWSSN  records  of  intermediate  focus  earthquakes  in  the 
Tonga,  Kermadec,  and  Kuril  regions  and  of  nuclear 
explosions  and  shallow  earthquakes  in  the  Aleutian 
region.  These  amplitudes  were  corrected  for  source 
mechanism.  The  Aleutian  data  were  sufficient  to  show 
that  intermediate  focus  earthquakes  in  that  region  occur 
in  the  colder  regions  of  the  slab.  At  short  periods 
shadowing  effects  which  could  be  associated  with  the 
slab  were  not  very  marked,  less  than  a  factor  of  2 
reduction  for  epicentral  distances  greater  than  50  degrees 
and  possibly  more  reduction  for  epicentral  distances 
between  30  and  50  degrees.  No  systematic  effects  due 
to  plates  were  found  in  the  long  period  data.  Some 
stations  in  the  predicted  shadow  zone  of  a  Tonga  earth¬ 
quake  recorded  low  amplitude  precursors  which  probably 
were  greatly  defocused  waves  which  ran  the  full  length 
of  the  slab.  Simple  diffraction  is  incapable  of  ex¬ 
plaining  the  short  period  results. 
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Numerical  and  analytic  models  were  constructed  to 
examine  four  hypotheses  for  the  origin  of  island  arc 
volcanics:  (1)  melting  at  temperature  lowered  by  inclu¬ 

sion  of  subducted  material;  (2)  frictional  heating  re¬ 
lated  to  the  descent  of  the  slab;  (3)  upwelling  of  mater¬ 
ial  from  the  asthenosphere  into  the  lithosphere  of  the 
island  arc;  (4)  concentration  of  pre-existing  melt  in 
the  asthenosphere.  Presently  available  geochemical  data 
are  insufficient  to  resolve  whether  a  significant  portion 
of  subducted  material  erupts  on  island  arcs.  This  does 
not  create  a  difficulty  in  evaluating  the  other  hypotheses, 
as  the  observed  eruption  of  temperature  of  island  arc 
volcanics  is  similar  to  that  of  basalt.  The  viscosity  of 
the  melt  and  the  geometry  of  the  asthenosphere  are  not 
conducive  to  segregation  of  melt.  At  low  concentrations 
of  melt,  the  regions  with  the  highest  concentration  of 
melt  contribute  disproportionately  to  the  magma  which 
segregates.  About  4  kb  of  shear  stress  is  needed  to 
cause  melting  above  the  slab.  Frictional  heating  is  likely 
to  cause  the  widespread  high  crustal  temperature  observed 
beneath  island  arcs.  Hypothesis  3  involves  no  obvious 
thermal  or  mechanical  difficulties  although  the  uncertain 
rheological  properties  of  the  region  near  the  base  of 
the  lithosphere  preclude  the  computation  of  a  realistic, 
detailed  model.  However,  it  is  both  mechanically  and 
thermally  reasonable  that  the  slab  could  entrain  inter¬ 
mediate  viscosity  material  near  the  base  of  the  lithosphere. 
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An  influx  of  material  from  the  asthenosphere  to  replace 
the  entrained  material  could  produce  high  crustal  temper¬ 
atures  and  provide  a  source  region  for  island  arc 
volcanics . 

The  hypotheses  and  models  studied  in  relation  to 
island  arc  volcanism  are  also  relevant  to  spreading  in 
intra-arc  basins.  Frictional  heating  above  the  slab 
is  probably  insufficient  to  cause  this  spreading.  A 
numerical  study  was  made  to  see  if  the  extension  could 
be  related  to  viscous  flow  of  material  entrained  with 
the  slab.  This  calculated  flow  had  the  correct  geometry 
to  cause  tension  behind  the  island  arc.  Intra-arc 
spreading  would  be  most  likely  to  occur  if  the  flow 
was  induced  in  the  intermediate  viscosity  region  be¬ 
tween  the  lithosphere  and  the  asthenosphere. 

3 . 8  Evolution  of  the  Downgoing  Lithosphere  and  the 

Mechanisms  of  Deep-Focus  Earthquakes  by  M.N.  Toksttz, 
N.H.  Sleep  and  A.T.  Smith 
Summary 

The  thermal  evolution  of  the  lithospheric  slab  at 
subduction  zones  and  its  geophysical  effects  are  numeri¬ 
cally  calculated.  An  alternating-direction,  implicit, 
finite-difference  scheme  is  used  to  compute  the  thermal 
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models  taking  into  account  all  heating  sources  and  phase 
boundaries.  These  models,  with  the  appropriate  spreading 
rates  and  dip  angles,  are  compared  with  different  island 
arc  systems.  Temperatures  inside  the  slab  are  strongly 
controlled  by  the  conductivity  and  by  the  time  elapsed 
since  the  initiation  of  descent.  The  depth  to  which 
temperature  anomalies  persist  is  generally  about  700  km  or 
less . 

The  thermal  results  are  used  to  construct  the  seismic 
velocities  and  ray  paths,  density  anomalies,  and  the  re¬ 
sulting  stress  distribution.  Comparing  the  theoretical 
stress  distribution  and  the  focal  mechanism  studies  of 
intermediate-  and  deep-focus  earthquakes  indicates  the 
importance  of  both  the  mantle's  rheology  and  the  tempera¬ 
ture-dependence  of  the  slae's  elastic  properties.  The 
intermediate  and  deep  focus  earthquakes  are  located 
along  the  coolest  region  of  the  slab.  The  theoretical 
results  explain  the  source  mechanisms  and  the  orientation 
of  principal  stresses  under  major  island  arcs. 
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INTRODUCTION 

The  interest  in  the  distribution  and  nature  of  deep- 
focus  earthquakes  has  been  intensive  and  the  circum-Pacif ic 
belt  has  been  extensively  studied  (Honda  1934;  Wadati  1935; 
Gutenberg  &  Richter  1939,  1954;  Oliver  &  Isacxs  1967;  Sykes 
1966;  Kondorskaya  &  Postolenko  1959).  A  renewed  interest  in 
and  much  better  understanding  of  island  arcs  and  deep 
seismicity  have  emerged  in  recent  years  with  the  sea-floor 
spreading  and  the  global  tectonic  hypotheses  (Isacks  et  al. 
1971;  Griggs  1972;  Smith  &  Toksttz  1972) . 

In  this  paper  we  discuss  the  thermal  regime  and  the 
stress  field  inside  a  descending  lithospheric  plate.  First 
we  present  theoretical  models  of  the  thermal  evolution  of  a 
downgoing  slab.  Later,  we  c  ompute  the  resulting  stress  fields 
and  compare  these  with  the  distribution  and  mechanisms  of 
intermediate  and  deep  focus  earthquakes. 
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THERMAL  EVOLUTION  OF  A  DOWNGOING  SLAB 

The  motion  of  the  descending  slab  is  a  continuation  of 
the  lithospheric  spreading.  Once  the  kinematics  of  the  motion 
are  specified,  the  calculation  of  the  temperature  field 
becomes  practical  without  solving  the  full  convection  problem 
for  the  earth. 

Studies  of  the  thermal  regimes  of  downgoing  slabs  have 
been  carried  out  by  several  investigators  using  different 
models  and  techniques.  McKenzie  (1969)  used  a  two-dimensional 
analytic  formulation  which  assumed  the  top  and  the  bottom  of 
the  slab  were  isotherms.  Griggs  (1972)  included  the  phase 
changes  in  a  numerical  solution  with  similar  boundary  con¬ 
ditions.  Turcotte  &  Oxburgh  (1969)  obtained  a  solution  using 
boundary  layer  theory.  Minear  &  Toksflz  (1970a, b),  Hasebe  et  al. 
(1970),  and  Toksttz  et  al.  (1971),  used  two- 

dimensional  numerical  models  in  which  the  mantle  surrounding 
the  slab  was  explicitly  included  so  energy  was  conserved. 

In  this  section  we  briefly  review  the  numerical  calcu¬ 
lations  following  Toksftz  et  al.  (1971)  . 

Wo  especially  focus  our  attention  to  effects  of  various 
physical  properties  and  to  the  questions  related  to  geometry 
and  the  rate  of  descent. 
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Computation  of  Temperature  Field 

The  computational  scheme  used  in  this  study  is  similar 
to  the  method  described  in  detail  in  an  earlier  paper  (Minear 
&  Toksttz  1970a) .  The  basic  model  consists  of  a  slab  of 
material  moving  downward  into  the  mantle  at  a  specified  angle. 
The  surrounding  mantle  material  is  assumed  to  be  fixed  (Fig. 

1)  .  The  computational  scheme  consists  of  translating 
temperatures  downward  in  the  slab  and  then  allowing  the  slab 
to  warm  up  over  the  time  interval  At,  which  corresponds  to 
the  vertical  movement  d.  In  essence  we  assume  the  dynamics 
and  compute  the  temperature  field  given  this  motion  field. 

Temperatures  are  computed  from  the  conservation  of 
energy  equation 

3T 

Vat  "  v  ’  <KVT>  +  h  at 

where 

Cp  =  specific  heat  at  constant  pressure,  p  =  density, 

T  =  temperature,  K  =  conductivity,  and  H  =  heat 
generation  rate  per  unit  volume. 

The  alternating-direction,  implicit,  finite-difference  scheme 
(Peaceman  &  Rachford  1955)  is  used  to  solve  (1)  numerically. 
The  method  was  described  earlier  by  Minear  &  Tokstiz  (1970a). 


Physical  Parameters  and  Energy  Sources 


Certain  physical  parameters  (conductivity,  specific 
heat,  density,  mantle  geotherm)  and  the  energy  sources  that 
are  incorporated  in  H  (radioactivity,  adiabatic  compression, 
phase  changes  and  shear-strain  heating)  must  be  specified 
for  solving  (1)  .  Some  parameters  can  be  specified  and 
remain  nearly  constant  as  a  function  of  position.  Others 
strongly  depend  on  temperature  and  need  to  be  updated  at 
each  step  and  every  grid  point. 

Physical  Parameters 

In  these  calculations  we  used  the  average  density  profile 
of  the  oceanic  model  (Press  1970)  and  a  constant  specific  heat 
of  Cp  =  1.3  x  107  erg/g  °C . 

In  the  thermal  conductivity  the  contributions  of  both 
the  lattice  conduction  and  radiative  heat  transfer  were  taken 
into  account.  One  set  of  models  was  based  on  the  formulation 
of  MacDonald  (1959)  . 

k  =  +  (16n3sT3)/(neQ  +  1207ro0exp  (-E/KT) )  (2) 

where 

k!  =  lattice  conductivity,  0.25  x  106  erg/cm-sec- °C , 

T  =  absolute  temperature,  K  =  Boltzmann's  constant, 
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E  =  width  of  energy  gap  for  electronic  conduction,  3  e.v., 
s  =  Stefan-Boltzmann  constant,  n  =  index  of  refraction 
(n  =  1.7)  ,  Gq  =  low  temperature  opacity  Kg  =  10  cm  ^)  , 
and  0^  =  electrical  conductivity  (Oq  =  10  ohm  ^cm  . 


The  other  set  of  models  used  the  experimental  results  of 
Schatz  &  Simmons  (1972)  on  olivine.  In  this  case  the 
lattice  conductivity  is  the  larger  of 


k1  =  (30.6  +  0.21T) 
k±  =  0.003  +  (3  x  10_6)z 


(3) 


where  z  is  in  kilometers,  T  in  °K  and  the  unit  of  conduc¬ 
tivity  is  cal/cm-sec-°C .  The  radiative  conductivity  is  given 
by 

kr  =  0  T  -  500°K 

(4) 

kr  =  5.5  x  10“6(T  -  500)  T  >  500°K 

The  conductivity  change  associated  with  the  phase  change 
at  400  km  is  probably  not  very  large  since  ferrous  iron, 
which  increases  opacity  and  reduces  radiative  conductivity, 
erters  the  common  phases  on  both  sides  of  the  transition  zone 
(Ringwood  1970) .  If  iron-free  oxides  such  as  MgO,  Al20^ 
or  Si02  exist  as  minerals  below  the  600  km  discontinuity, 
the  radiative  conductivity  would  be  extremely  high. 

In  earlier  studies  the  unperturbed  geotherm  of  the  mantle 
was  either  assumed  to  be  adiabatic  (McKenzie  1969) ,  the 
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conductivity  geotherin  (Minear  &  Toksflz  1970a, b;  Tokstiz  et  al. 
1971)  or  to  be  controlled  by  the  solidus  temperature  of 
peridoti te  (Griggs  1972) .  In  this  paper  an  approximation  to 
the  average  geotherm  in  a  convecting  earth  was  also 
used.  The  unperturbed  geotherm  is  indicated  in  each  tempera¬ 
ture  figure  in  the  following  section. 

Heat  Sources 

In  addition  to  conduction  of  heat  from  the  surrounding 
mantle,  the  slab  is  heated  by  internal  heat  sources.  These 
consist  of  radioactivity,  adiabatic  compression,  phase  changes, 
and  the  shear  strain  heating  along  the  slab-mantle  boundaries. 

Radioactive  neating .  The  average  radioactivity  of  the  oceanic 
lithosphere  is  probably  about  equal  to  the  radioactivity  of 
the  upper  mantle,  because  differentiation  at  mid-oceanic 
ridges  occurs  mainly  above  30  km  depth  (Kay  et  al. 1970) 

Radioactive  heating  in  the  upper  mantle  of  the  earth  is 
too  low  to  seriously  affect  the  evolution  of  the  slab.  A 
more  important  effect  of  radioactive  heating  is  on  the  un¬ 
perturbed  geotherm  of  the  upper  mantle.  The  estimates  of 
heat  generation  due  to  radioactivity  cover  a  wide  range: 

_7 

2.3  x  10  erg/gm  sec  (MacDonald  1959,  1963),  5.4  x  10~8 
erg/gm  sec  (Armstrong  19G8) ,  and  1.5  x  10_8  erg/gm  sec 
(Hurley  1968a, b).  The  larger  radioactivity  of  MacDonald  would 
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increase  the  temperature  only  6°C  in  the  10  m.y.  descent  of 
the  slab.  In  these  calculations  we  have  used  heat  pro¬ 
duction  values  of  2.29  x  10  ^  ergs/g  sec  for  the  upper 

“  8 

mantle  and  1.67  x  10  ergs/g  sec  for  the  mantle  below  465  km. 


Adiabatic  heating.  As  the  lithosphere  descends  into  the  mantle, 
it  is  compressed  and  heated.  The  adiabatic  temperature  gradient 
is  given  by 


3T/3z  =  gaT/Cp  (5) 

where  g  is  the  gravitational  acceleration  and  a  is  the 
volume  coefficient  of  thermal  expansion.  The  rate  of  thermal 
energy  release  at  depth  h  due  to  adiabatic  compression  is 


dQ 

dt 


h 


=  CpP(3T/3t) 


=  C  p (3T/3z) V  =  pgaTV 
p  z  z 


(6) 


where  Vz  is  the  vertical  component  of  the  velocity  of  the 
downgoing  slab.  The  value  dQ/dt  can  be  evaluated  at  each 
point  by  using  appropriate  values  for  p,  g,  a,  and  T  (Hanks 
&  Whitcomb  1971) .  p  and  g  are  known  and  T  is  computed  from 
the  previous  time  step.  Laboratory  measurements  of  a  tabu¬ 
lated  by  Skinner  (1966)  indicated  a  slight  increase  in  a  with 
increasing  temperature  and  a  decrease  in  a  with  increasing 
pressure.  For  these  calculations  we  expressed  the  depth 


dependence  of  a  by 

a  =  exp(3.58  -  0.0072z)  (7) 

where  z  is  depth  in  km  and  a  is  measured  in  10~6  °C_1.  Near 
the  surface  (z  =  0)  the  value  of  ot  =  36  x  10  ^  °C  ^  corresponds 
to  that  of  olivine  at  about  T  =  500°C  (Skinner  1966) .  The 
exponential  coefficient  0.0072,  for  depth  dependence,  is 
somewhat  lower  than  estimates  of  Birch  (1968)  and  slightly 
higher  than  those  of  Verhoogen  (1951) .  The  above  expression 
gives  a  value  of  a  =  18  x  lO-6  °C_1  at  z  =  1000  km,  in 
agreement  with  the  value  used  by  Griggs  (1972)  . 

Since  a  decreases  with  depth  while  the  temperature  is 
increasing  inside  the  slab,  the  adiabatic  heat  generation 
(dQ/dt)  varies  slowly,  increasing  slightly  with  depth. 

Phase  changes.  Several  phase  changes  may  occur  in  the  upper 
mantle.  Seismic  velocity  profiles  indicate  the  presence  of 
second-order  discontinuities  at  depths  of  about  350  and  650 
km  (Toksttz  et  al.  1967  ;  Johnson  1967;  Julian  &  Anderson  1968; 
Archambeau  et  al .  1969)  that  are  most  likely  caused  by 
olivine  to  spinel  and  spinel  to  post-spinel  changes  (Anderson 
J967).  In  addition  to  these,  a  phase  change  corresponding  to 
a  basalt-eclogi te  or  a  plagioclase  peridoti te-garnet  peridotite 
teaction  is  considered  at  shallower  depth. 

If  the  vertical  velocity  through  a  phase  change  is  fast 


i. 
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enough  to  disturb  the  isotherms,  the  phase  boundary  will  be 
moved  from  its  normal  position  (Schubert  &  Turcotte  1971) . 

To  calculate  the  effects  of  phase  changes, all  points  in  the 
model  were  examined  at  the  end  of  each  translation  step  and 
the  change  in  the  amount  of  each  phase  from  the  previous  step 
recorded.  Latent  heat  was  calculated  from  this  change  and 
included  as  a  heat  source  in  the  model.  The 

routine  for  translating  temperature  was  also  used  to  determine 
the  amount  of  each  phase  with  the  slab.  This  calculation 
scheme  was  used  for  phase  changes  going  both  ways  whether  or 
not  the  slab  is  assumed  to  move. 

Three  phase  changes  were  considered  in  the  models. 

Phase  change  energies  were  computed  from  entropy  (S)  and 
volume  (V)  changes.  Because  the  temperatures  vary  inside 
the  slab,  the  phase  boundaries  were  determined  from  the  dP/dT 
slopes  for  each  phase  change,  starting  from  the  shallowest. 

The  AS  values  for  the  three  phase  changes  are  taken  as 
-11.7  x  10  ,  -7.13  x  10  ,  and  -5.94  x  10^  ergs/g°C.  Corre¬ 
sponding  AV  values  are  -6.5  x  lO-2,  -2.65  x  10-2,  and  -2.62  x 
-2  3 

10  cm  /g.  These  values  are  averages  of  those  given  by 
Vorhoogen  (1965) ,  Akimoto  &  Fujisawa  (1968) ,  Ringwood  (1970, 
1972),  Solar  et.  a_L.  (  1964).  The  parameters  for  the  spinel- 
post-spinel  phase  changes  are  the  least  reliable.  These 
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reactions  were  studied  indirectly  with  chemical  analogs 
(Akimoto  1970;  Ringwood  &  Major  1970) .  Ringwood  (1972) 
takes  this  reaction  to  be  exothermic  (and  dP/dT  >  0)  while 
others  (D.L.  Anderson  1972,  personal  communication)  suggest 
that  it  may  be  endothermic.  Because  of  these  uncertainties, 
we  constrained  the  phase  boundary  to  remain  at  650  km  depth, 
but  adopted  a  positive  heat  contribution.  We  also  calcu¬ 
lated  some  models  without  any  heat  contribution  from  this 
phase  boundary.  No  allowance  was  made  for  any  changes  in 
heat  balance  due  to  possible  melting  of  the  material. 

Shear-strain  heating.  An  accurate  estimate  of  the  shear- 
strain  heating  is  extremely  difficult,  because  we  know  neither 
the  mechanisms  involved  in  the  descent  of  the  slab  nor  the 
parameters,  such  as  the  thickness  of  the-  shear  zone  and  its 
viscosity.  The  complexity  of  the  problem  even  for  an  ideal¬ 
ized  fluid  model  has  been  illustrated  by  Turcotte  &  Oxburgh 
(1969b)  . 

In  the  calculations  we  used  several  values  for  shear- 
strain  heating.  In  all  cases  viscous  heat  generation  was 
confined  to  14  km  thick  layers  along  the  top  and  bottom  of 
the  slab.  At  the  top  edge,  the  maximum  shear-strain  heat 
generation  rates  are  1.6  x  10  ergs/cm^  sec,  unless  specified 
otherwise.  At  the  lower  boundary,  the  shear-strain  heating 

3 

ergs/cm  sec  (a  factor  of  ten  less  than  the  top) 


is  1.6  x  10 


and  it  has  a  negligible  effect  on  the  thermal  regime. 


At  the  8  cm/yr  subduction  rate,  the  frictional  heating 
used  at  the  top  of  the  slab  is  equivalent  to  a  stress  of 
a  few  kilobars  (in  an  extreme  case  it  could  be  as  high  as 
4  kb)  .  Although  this  is  high  compared  to  apparent  stress 
drop  associated  with  earthquakes,  there  is  evidence  that 
high  stresses  prevail  in  the  lithosphere  at  the  convergence 
zones.  For  example,  to  produce  the  observed  uplift  at  the 
outer  rise  of  the  Japan  trench,  stresses  of  a  few  kilobars 
are  needed  (Hanks  1971)  . 

Temperature  Models 

A  series  of  temperature  models  were  computed  to 
demonstrate  the  evolution  of  the  thermal  regime.  Only  a 
few  models  will  be  shown  in  this  paper  in  addition  to  those 
in  the  previous  papers  to  demonstrate  the  effects  of 
important  parameters  and  varying  geometry  (Minear  &  Toksttz 
1970a, b;  Toksttz  et  al.  1971) . 

The  development  of  the  thermal  regime  in  a  downgoing 
slab  is  shown  in  Figs.  2  and  3  as  a  function  of  time.  Temperature 
fields  shown  after  3.6  and  7.1  m.y.  (Fig.  2)  and  10.7  m.y. 

(Fig.  3)  from  the  start  of  downward  motion  demonstrate  the 
relative  effects  of  time  and  crossing  of  phase  boundaries. 

The  phase  boundaries  are  elevated  by  as  much  as  150  km  as  a 
result  of  the  lower  temperatures.  At  the  rate  of  8  cm/yr,  the 
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interior  of  the  slab  remains  cooler  relative  to  the  surroundings 
to  a  depth  of  at  least  600  km.  In  the  coolest  zone,  the 
temperature  difference  is  more  than  500°C. 

The  rapid  heating  of  the  slab  interior  below  600  km 
is  primarily  due  to  the  higher  thermal  conductivity  (because 
of  increased  contribution  of  radiative  heat  transfer  at 
higher  temperatures)  and  to  the  contribution  of  spinel- 
post-spinel  phase  transformation.  The  thermodynamic  param¬ 
eters  of  this  transformation  are  not  well  known  and  that  is 
why  we  constrained  the  phase  boundary  to  a  constant  depth. 
However,  if  the  reaction  is  exothermic,  the  heat  generated 
tends  to  raise  the  temperature  to  the  ambient.  When  the 
internal  slab  temperature  approaches  the  mantle  temperature, 
the  heat  generated  internally  cannot  be  transferred  to  cooler 
portions  and  the  source  region  temperature  increases  rapidly. 

As  a  result  the  slab  loses  its  integrity  and  is  assimilated 
into  the  mantle.  Even  without  any  heating  from  the  lower 
phase  boundary,  the  slab  reaches  thermal  equilibrium  at  about 
750-800  km  depth  (Toksflz  et  al.  1971). 

Parameters  that  strongly  influence  the  temperatures 
inside  and  around  the  downgoing  lithosphere  include  the 
conductivity,  and  the  descent  rate.  The  effect  of  using  the 
Schatz  &  Simmons  (1972)  conductivity  is  shown  in  Fig.  4. 

Comparing  this  with  Fig.  3  [based 


on  conductivity  given  by 
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(2)]  illustrates  that  in  the  case  of  lower  effective 
conductivity  (Fig.  4)  temperatures  are  lower  inside  the  slab. 

The  temperatures  inside  the  slab  are  strongly  controlled 
by  the  time  elapsed  after  the  initiation  of  the  descent.  A 
slower  moving  slab  will  tend  to  heat  up  more  because  of 
the  increased  conduction  of  heat  from  the  surrounding  mantle. 
For  a  descent  velocity  corresponding  to  1  cm/yr  spreading 
rate  and  a  total  time  of  100  m.y.,  Toksftz  et  a^.  (1971) 

showed  that  with  parameters  similar  to  those  of  Fig.  3  the 
slab  reaches  thermal  equilibrium  at  a  depth  of  about  400  km. 
Here  we  expanded  these  calculations  to  include  stop-and-start 
type  slab  motions  and  differing  angles  at  which  the  slab 
penetrates  into  the  mantle. 

The  temperature  field  of  a  slab  dipping  29  degrees  was 
computed  using  the  physical  parameters  of  the  base  model 
(i.e.,  Fig.  3)  and  the  results  are  shown  in  Fig.  5.  Although 
16  m.y.,  as  opposed  to  11  m.y.  for  the  base  model,  was  re¬ 
quired  for  this  slab  to  penetrate  to  650  km,  the  maximum 
depth  of  penetration  of  both  slabs  is  controlled  by  the  lower 
phase  change.  The  maximum  penetration  of  mantle  isotherms 
into  the  slab  is  about  30  to  100  km  less  in  the  shallower 
dipping  model. 


It  is  conceivable  that  for  a  period  of  time  a  slab  may 
become  hung  up  at  a  trench  such  that  no  subduction  occurs. 
Once  subduction  ceases  it  is  likely  that  the  slab  will 
detach  and  sink  as  it  may  have  done  in  che  New  Hebrides. 
Detachment  could  also  occur  if  new  material  were  being 
subducted  more  slowly  than  the  separation  rate.  The  thermal 
effects  of  stopping  the  motion  after  the  slab  penetrated 
to  650  km  are  shown  in  Fig.  6.  In  these  calculations  the 
same  slab  shown  in  Fig.  5  was  allowed  to  remain  stationary. 
After  16  m.y.  (total  age  32  m.y.)  the  slab  was  still  evident 
but  the  temperature  anomaly  was  spread  over  a  larger  region. 
After  48  m.y.  only  a  very  broad  temperature  anomaly  remained. 
It  is  not  likely  that  a  seismic  zone  could  remain  active 
for  a  long  time  after  subduction  ceases.  The  most  likely 
event  is  that,  once  the  slab  is  partially  heated,  detachment 
would  occur  and  the  cooler  slab  would  sink. 

The  use  of  a  cooler  ("convective")  geotnerm  and  the 
effect  of  the  600  km  phase  boundary  is  shown  in  Fig.  7.  The 
use  of  the  "convective"  geotherm  does  not  greatly  alter  the 
temperature  field  (Fig.  7).  The  lower  phase  change  was  not 
included  in  constructing  the  gootherm  of  the  model.  Unlike 
the  models  using  the  strongly  superadiabatic  MacDonald  (1959) 
geotherm,  the  temperature  fieM  is  not  greatly  influenced 
by  the  kinematics  assumed  for  the  slab,  since  motions  in  the 
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nearly  adiabatic  regions  of  the  mantle  effect  temperatures 
only  slightly.  It  should  also  be  noted  that  the  effect  of  a 
400  km  phase  change  is  evident,  although  this  reaction  was 
explicitly  included  when  constructing  the  geotherm.  The  600 
km  phase  change  was  not  included  in  these  calculations. 

The  omission  of  the  phase  change  at  600  km  depth  in¬ 
creases  the  depth  at  which  the  slab  reaches  thermal 
equilibrium.  Even  in  this  case,  however,  the  equilibration 
occurs  at  about  800  km  depth  (Toksttz  et  al .  1971) .  At  these 
great  depths,  the  contribution  of  radiative  heat  transfer 
becomes  substantial  and  the  slab  interior  warms  up  rapidly. 
Continuing  the  calculations  to  greater  depths  results  in 
temperatures  inside  or  around  the  slab  higher  than  the  ambient 
mantle  temperatures.  It  is  reasonable  to  assume  that  the 
slab  assimilates  into  the  mantle  and  becomes  part  of  the 
general  mantle  convection  pattern  after  it  reaches  thermal 
equilibrium.  As  long  as  the  subduction  continues  at  some 
uniform  rate,  the  temperatures  inside  the  slab  above  the  depth 
at  which  the  slab's  interior  reaches  equilibrium  change 
very  little  with  time.  For  an  8  cm/yr  subduction  rate  for 
example,  the  thermal  regime  reaches  a  nearly  stable  state 
alter  about  10  or  12  million  years  above  a  depth  of  700  km. 

divergence  of  continental  lithospheres.  In  a  separate  yet 
related  area  to  the  subduction  of  the  oceanic  lithosphere,  we 
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mvestigated  the  thermal  regime  for  the  convergence  of  two 
continental  lithospheres.  Because  of  the  generally  shallower 
angle  of  penetration  [for  example,  the  dip  of  the  thrust 
belt  is  between  10  and  20°  for  the  Zagros  convergence  zone 
in  Southern  Iran  (Nowroozi  1972)],  slower  rate  of  subduction, 
continental  geotherm,  and  higher  radioactivity  in  the 
continental  crust,  the  thermal  evolution  is  expected  to  be 
different  from  that  in  the  oceanic  case. 

The  large  scale  structure  in  the  region  of  a  continent- 
continent  collision  may  be  in  part  inherited  from  the  pre¬ 
existing  active  continental  margin  and  in  part  due  to  the 
subducted  continental  crust.  For  continental  crust  to  have 
much  effect  on  temperature  during  its  subduction,  either  the 
shear  strain  heating  would  have  to  be  high  or  the  subduction 
rate  would  have  to  be  low  for  prolonged  radioactive  heating. 
The  average  radioactive  heat  production  of  continental 
crust  is  about  4.1  x  10  6  erg/gm-sec  (Hurley  1968a, b; 
Armstrong  1968)  and  this  can  cause  a  maximum  temperature 
increase  of  10°C/m.y. 

A  numerical  model  of  continental  convergence  is  shown 

in  Fig.  8.  The  subduction  rate  is  1  cm/yr.  A  high  value 

-4  3 

of  frictional  heating  (1.6  x  10  ergs/cm  sec)  is  used 

along  the  upper  edge.  Radioactive  heat  generation  is 

4.1  x  10  6  erg/gm-sec  in  the  crust  (hatched  region)  and 
—  8 

1.5  x  10  erg/gm-sec  in  the  mantle.  Although  this  model 


shows  high  geothermal  gradients  above  a  subducted  continent 
this  does  not  seem  to  be  sufficient  to  cause  extensive 
melting  or  magmatization  associated  with  orogenies. 

Another  continental  convergence  model  with  lower 
frictional  heatings  (1  x  10  ^  ergs/cm^  sec  above  35  km  depth) 
is  shown  in  Fig.  9  for  two  time  intervals.  Other  parameters 
are  the  same  as  those  of  Fig.  8.  After  7  7  m.y.  the 
temperatures  are  slightly  lower  inside  the  slab,  except  at 
the  boundary  below  100  km  depth.  After  about  50  m.y. 

(lower  figure)  the  temperature  returns  to  steady  state. 

We  may  conclude  that  orogeny  is  most  likely  to  occur  during 
the  collision  or  after  about  50  m.y.  when  the  subducted 
continental  material  has  become  heated.  In  the  former  case 
the  heat  of  the  orogeny  would  be  related  to  volcanic  activity 
along  the  active  margin  prior  to  the  collision.  In  the 
latter  it  would  be  radioactive  decay  in  the  subducted  material. 


Geophysical  Effects 

The  strongly  varying  temperature  field  of  the  upper 
mantle  inside  and  in  the  vicinity  of  the  slab  affects  the 
surface  heat  flux,  volcanism,  density,  seismic  wave  velocity 
and  attenuation,  as  well  as  the  stress  field  and  the  occurrence 
and  mechanisms  of  earthquakes.  Most  of  these  effects  have 
been  discussed  in  previous  papers  (Toksttz  et  aj^.  1971;  Sleep 
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1973) .  Not  all  surface  data  are  equally  revealing  about  the 
properties  of  the  deep  lithosphere.  For  example,  heat  flow 
is  primarily  controlled  by  the  temperature  field  at  relatively 
shallow  depth  and  by  convective  heat  transfer  underneath  and 
behind  island  arcs  (Hasebe  et  al.  1970) . 

The  gravity  anomalies  associated  with  downgoing  slabs 
are  broad  regional  anomalies  (Minear  &  Toksttz  1970a, b; 

Griggs  1972) .  These  are  generally  obscured  by  the  much 
sharper  anomalies  over  trenches  and  island  arcs  that  result 
from  topography  and  crustal  structure.  Observed  gravity 
anomalies  have  been  fitted  with  a  model  including  a  slab 
in  the  Aleutians  (Grow  1972)  . 

The  most  pronounced  effects  of  the  temperature  fields 
of  descending  slabs  are  on  seismic  velocities,  Q  structure 
and  intermediate  and  deep  focus  earthquakes.  Velocity  and 
amplitude  anomalies  of  seismic  waves  have  been  used  to 
demonstrate  the  existence  of  the  slab  (Davies  &  McKenzie  1969; 
Sorrells  et  al.  1971;  Mitronovas  &  Isacks  1971;  Toksflz  et  al. 
1971;  Jacob  1970,  1972;  Abe  1972a, b;  Davies  &  Julian  1972). 

The  observed  travel  time  anomalies  are  in  close  agreement 
with  the  theoretical  models  based  on  the  calculated  temperature 
profiles  (Sleep  1973) . 

The  distribution  and  stress  fields  associated  with  these 
deep  focus  earthquakes  will  be  discussed  in  detail  in  the 


next  section. 
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DISTRIBUTION  AND  MECHANISMS  OF  DEEP  FOCUS  EARTHQUAKES 

Since  their  discovery  the  origin  and  mechanisms  of 
intermediate  and  deep  focus  earthquakes  have  been  of  great 
interest.  The  sea-floor  spreading  and  global  tectonic 
hypothesis  with  descending  slabs  has  provided  an  explanation 
for  these  deep  earthquakes.  It  is  generally  accepted  that 
the  slab  behaves  like  a  stress  guide.  However,  the  mechanisms 
of  the  intermediate  and  deep  earthquakes  show  remarkable 
differences  between  regions  not  explained  by  any  simplified 
hypothesis  (Isacks  &  Molnar  1971) .  In  this  section  we 
treat  this  problem  in  two  steps.  First  we  determine  the 
hypocenters  relative  to  the  slab  boundaries.  Then,  computing 
the  stresses  arising  from  subduction,  we  investigate  the 
mechanisms  of  these  earthquakes. 

Distribution  of  Deep  Focus  Earthquakes 

The  shallow  earthquakes  near  island  arcs  are  concen¬ 
trated  in  the  convergence  zones  between  plates,  mostly  along 
the  upper  boundary  of  the  underthrusting  plane  of  the  oceanic 
lithosphere.  For  the  Aleutians  and  Japan  both  the  distribution 
of  hypocenters  and  the  source  mechanisms  of  earthquakes 
support  this  conclusion  (Barrazangi  &  Dorman  1969;  Stauder 
1968;  Utsu  1971;  Isacks  &  Molnar  1971;  Plafker  1972;  Engdahl 
1972) .  for  intermediate  and  deep  focus  earthquakes  the 
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precise  location  of  hypocenters  becomes  more  difficult 
because  of  the  effect  of  velocity  anomalies  within  the  slab 
on  the  seismic  ray  paths. 

to  further  constrain  the  locations  of  inter¬ 
mediate  earthquakes  with  respect  to  the  slab,  we  investigated 
the  central  Aleutians.  Nuclear  explosions  LONGSHOT,  MILROW 
and  CANNIKIN  provided  sources  of  precise  location  and  time  for 
calibrating  the  velocity  models  and  slab  location  relative  to 
the  trench.  Furthermore  the  network  of  seismic  stations  in 
the  area  provided  very  good  locations  of  the  earthquakes. 

The  theoretical  ray  paths  from  surface  focus  events  were 
calculated  using  seismic  velocities  deduced  from  a  temperature 
model  of  the  A]  utian  slab.  A  most  likely  location  of 
LONGSHOT  relatx  e  to  the  slab  was  determined  by  comparing 
the  computed  travel  time  delays  (Fig.  10)  and  shadow  zone 
with  the  observations  (Abe  1972a;  Jacob  1972;  Sleep  1973). 

The  preferred  slab  location  based  on  the  above  data  is  then 
compared  with  the  hypocenters  of  earthquakes  (Engdahl  1971) 
in  Fig.  11.  It  is  very  clear  that  the  intermediate  focus 
earthquakes  are  located  along  the  coolest  region  of  the  slab. 
Similar  distributions  are  implied  for  the  earthquakes  under 
Japan  (Katsumata  1967;  Utsu  1971)  and  the  Tonga  Kermadec 
region  from  travel  times  and  attenuation  characteristics  of 
P  and  S  waves  (Sykes  e_t  al .  1969;  Mitronovas  et  al .  1969; 

Toksfiz  et  al.  1971;  Sleep  1973). 
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Mechanisms  of  Deep  Focus  Earthquakes 

Observational  data  on  source  mechanisms  of  intermediate 
and  deep  focus  earthquakes  have  accumulated  in  recent  years 
from  fault-plane  studies.  A  comprehensive  summary  is  given 
by  Isacks  &  Molnar  (1971) .  In  this  section  we  describe  the 
calculation  of  stress  field  in  a  downgoing  slab  and  compare 
this  with  the  earthquake  distribution  and  mechanisms. 

Computation  of  stress  fields.  The  theoretical  studies  of 
stress  regime  (McKenzie  1969;  Griggs  1972;  Smith  &  Toksttz 
1972)  are  hampered  by  both  the  computational  difficulties 
and  the  limited  knowledge  of  mantle  rheology.  Here,  following 
Smith  &  Toksttz  (1972)  we  utilized  a  kinematic  approach  for 
computing  stresses  in  a  subducting  slab.  We  included  the 
direct  effects  of  subduction  (i.e.,  body  forces  due  to 
density  anomalies  resulting  from  cooler  temperatures)  in  the 
calculations . 

In  our  calculations  a  slice  of  mantle,  which  includes  the 
slab,  was  modelled  as  a  heterogeneous,  temperature-dependent, 
viscoelastic  medium.  Using  the  correspondence  principle  for 
viscoelasticity  (Christensen  1971)  ,  the  operators  of  the  visco¬ 
elastic  problem  arc  replaced  by  an  equivalent  elastic  modulus 
which  incorporates  the  temperature  and  stress  dependence  for 
that  increment  in  time.  The  much  simpler  elastic  problem  was 
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then  solved  numerically  for  the  stresses  (Zienkiewicz  1971). 
For  sufficiently  small  time  increments,  this  also  applies  to 
a  nonlinear  problem.  The  model  for  the  slab  assumed  quasi¬ 
static  behavior  after  application  of  the  loads.  The  time 
since  loading  was  taken  as  approximately  100  years;  this 
approximates  the  frequency  of  transient  phenomena  (i.e., 
deep  earthquakes)  within  the  slab.  Thus  the  problem  assumes 
a  certain  temperature  dependent  rheology  and  application  time 
for  the  loads;  the  resulting  operator  E,  an  effective  Young's 
modulus,  is  solved  as  an  elastostatic  problem  for  the  strains 
resulting  from  creep  since  application  of  the  load. 

Using  the  thermal  models  in  the  preceding  section  and 
the  thermal  expansion  coefficient  (7) ,  the  density  anomalies 
were  computed  relative  to  the  mantle  isotherm  at  each  depth. 
These  specify  the  vertical  body  forces  (F  =  gpaAT)  at  all 
points  inside  the  slab.  For  a  temperature  model  similar 
to  Fig.  4,  the  corresponding  body  forces  are  shown  in  Fig.  12. 
The  density  changes  and  the  resulting  forces  due  to  phase 
changes  were  separately  taken  into  account. 

A  finite  difference  scheme  incorporating  the  body  forces 
was  then  used  to  compute  the  stresses  for  the  elastostatic 
analogy.  The  scheme  uses  an  integral  formulation  of  the 
equations  together  with  over-relaxation  to  obtain  a  solution 
(Tillman  1971;  Smith  &  Toksttz  1972);  thus,  it  is  essentially 
a  finite-element  computation.  A  20  km  grid  spacing  was  used 
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to  solve  the  problem  for  an  840  by  680  km  region  enclosing 
the  slab  and  the  adjacent  mantle.  The  imposed  boundary 
conditions  were  that  the  top  face  (the  earth  surface)  was 
free,  the  vertical  boundaries  at  the  sides  of  the  region 
and  the  bottom  boundary  were  specified  to  be  rigid.  Other 
boundary  conditions  were  also  used  for  the  side  boundaries, 
but  the  effects  were  found  to  be  insignificant  for  the 
calculated  stresses  inside  the  slab. 

The  rheology  of  the  slab  and  the  surrounding  mantle 
were  specified  in  terms  of  temperature  and  pressure  dependent 
constants.  For  the  mantle,  the  operator  E  was  obtained  using 
effective  mantle  viscosities  (McConnell  1968;  Cathles  1971) 
and  the  elastic-viscoelastic  analogy  (Christensen  1971) .  The 
slab  rheology  was  assumed  to  b  ;  temperature  dependent.  For 
long-term  stresses,  the  temperature  dependence  of  the  slab 
properties  are  most  likely  determined  by  creep  processes 
(Weertman  1970)  .  Two  separate  relationships  were  used  to 

express  the  variation  in  the  effective  Young’s  modulus  inside 
the  slab: 

Eslab  =  Emantle  +  5,0  x  1()1°  exp(-0.3S  x  10_2AT)  (8) 

and 

Eslab  -  Emantle  exPt-(1.2  x  10  2  +  5.5  x  10  5p)AT]  (9) 
where  AT  i.;  the  lateral  temperature  difference  between  the 
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mantle  and  slab,  and  p  is  the  pressure  in  kilobars.  The 
latter  represents  a  much  stronger  temperature  dependence. 
Numerical  models  calculated  with  both  of  these  relationships 
did  not  alter  the  stresses  significantly  (Smith  &  Toksttz 
1972) . 

The  marked  effect  that  the  mantle's  support  has  upon 
the  stress  distribution  inside  the  slab  is  illustrated  in 
Fig.  13.  Both  slab  models  are  identical  and  include  only 
the  body  forces  due  to  density  anomalies  resulting  from  the 
thermal  contraction.  The  stress  regimes  are  quite  different. 
With  little  support  or  resistance  to  the  descending  litho¬ 
sphere,  the  tension  axis  follows  the  trend  of  the  slab  along 
its  full  length.  On  the  other  hand,  when  the  mantle's 
support  increases  with  depth  (top  diagram.  Fig.  13),  the 
maximum  principle  stress  becomes  compressional  along  the  full 
length  of  the  slab. 

As  these  two  contrasting  models  indicate,  the  mantle 
support  and  viscosity  variation  with  depth  significantly 
affect  the  direction  of  stresses  and  the  foca]  mechanisms  of 
earthquakes.  For  short  slabs  where  deepest  penetration  does 
not  extend  below  low-velocity  or  the  low  viscosity  zone, 
the  direction  of  least  compression  (i.e.,  "tension")  axis 
should  parallel  the  slab's  axis.  For  long  slabs  where 
penetration  to  more  resistant  mantle  takes  place,  the  stresses 
become  compressive. 


-168- 


The  density  anomalies  due  to  phase  changes  introduce 
additional  stresses  in  the  region  of  the  transformation. 

Above  the  phase  change,  the  stress  is  down-dip  tension,  and 
below,  it  is  down-dip  compression.  For  the  o.1  ivine-spinel 
phase  change,  this  can  add  about  200  bars  of  shear  stress 
(Smith  &  Toksttz  1972) .  The  effect  of  these  phase  changes 
may  be  important  for  variations  in  deep  seismicity.  If 
the  post-spinel  phase  change  at  600  km  rises  within  the  slab, 
additional  stress  would  occur  below  the  transformation  and 
help  account  for  the  sudden  rise  in  seismicity  at  about 
600  km  observed  in  the  Tonga-Kermadec  region. 

The  stresses  from  convection  on  the  slab  in  the  mantle 
can  be  incorporated  as  surface  traction  without  including 
body  forces.  These  are  then  superimposed  on  the  previous 
solutions.  A  uniform  50  bar  resistance  along  the*  upper  face 
of  the  slab  introduces  primarily  compressive  stresses  parallel 
to  the  axis  of  the  slab.  The  resulting  shear  stress  does 
not  exceed  100  bars  throughout  the  length  of  the  slab  (Smith 
&  Toksttz  1972).  These  models,  however,  do  not  properly 
simulate  the  shallow  stresses  at  the  island  arcs,  and  that 
was  not  their  intent.  For  shallow  stresses  the  boundary 
conditions  at  the  sides  of  the  region  and  the  fault  zone 
between  the  two  lithospheric  plates  are  important  and  require 
further  refinements.  Yet  as  formulated  they  are  adequate 
for  the  problem  at  hand:  the  intermediate  and  deep  stresses 
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within  the  slab.  Thus  there  is  no  inconsistency  between  the 
relatively  high  stresses  (1  kb  or  more)  used  for  shear  heating 
calculations  and  the  models  of  the  stress  distribution  within 
the  slab's  interior.  To  resolve  those  problems  associated 
with  the  bending  and  thrusting  of  the  descending  lithosphere, 
additional  models  are  necessary. 

Comparison  with  earthquake  mechanisms.  The  comparison  of 
the  calculated  stress  fields  with  the  distribution  and 
mechanisms  of  intermediate  and  deep  focus  earthquakes  can  now 
be  made  in  regions  where  the  seismicity  is  best  understood. 

For  this  we  chose  Japan  (Honshu)  and  the  Kermadec  region. 

In  Fig.  14  the  theoretical  stresses  and  the  earthquake 
hypocenters  given  by  Utsu  (1971)  for  the  Japanese  island  arc 
are  compared.  The  shallow  earthquakes  are  concentrated 
primarily  along  the  shear  zone  between  the  underthrusting 
oceanic  lithosphere  and  the  continental  lithosphere.  As  far 
as  it  can  be  determined  from  travel  times  and  attenuation 
characteristics,  intermediate  and  deep  earthquakes  are 
located  in  the  interior  of  the  slab  along  a  band  corresponding 
to  the  maximum  shear  stresses . 

The  correspondence  between  the  theoretical  and  observed 
directions  of  maximum  stresses  is  shown  in  Fig.  15.  in  the 
Honshu  case  the  model  dips  30°  and  has  mantle  properties  as 
.shown  in  Fig.  J3a.  The  rheology  of  the  slab  is  defined  by 
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equation  (9)  .  On  the  right  the  directions  of  maximum  prin¬ 
cipal  stress  determined  by  the  earthquake  focal  mechanisms 
(Isacks  &  Molnar  1971)  are  shown.  Both  the  observed  and 
the  theoretical  maximum  stresses  are  down-dip  compression, 
which  are  in  good  agieement  with  each  other.  In  the  case 
of  Kermadec ,  where  tie  slab  has  a  higher  dip  and  rheology 
specified  by  equation  ^8),  the  observed  principal  stresses 
are  down-dip  compression  for  deep  earthquakes  and  "tension" 
for  intermediate  earthquakes.  Again,  the  observations  are 
in  good  agreement  with  the  theoretical  model. 

The  calculated  maximum  shear  stresses  inside  the 
descending  slab  are  approximately  500  bars.  Estimates  of 
stress  drop  are  available  for  a  limited  number  of  deep  focus 
earthquakes  and  are  near  100  bars  (Wyss  &  Molnar  1972) . 

This  figure  is  reasonable  for  an  initial  (in  situ)  stress 
of  about  500  bars  assuming  that  a  relationship  similar  to 
those  of  shallow  earthquakes  exists  between  the  initial  stress 
and  stress  drop. 

"Faulting"  mechanisms  of  deep  earthquakes.  The  fault  plane 
solutions  obtained  for  deep  focus  earthquakes  indicate  that 
double-couple  type  source  functirns  fit  the  observed  radiation 
pattern  (Isacks  &  Molnar  1971;  Stauder  1968) .  The  650  km 
deep  Colombian  earthquake  studied  in  detail  by  Mendiguren 
(1^72;  epicenter  1.5°S,  72.6°W,  origin  time  July  31,  1970, 
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17h  08m  05.4s,  magnitude  7.1)  using  both  first  motion 
polarities  and  free  oscillations  verifies  a  double-couple 
source.  Any  volumetric  component  of  this  earthquake  source 
contributes  no  more  than  10  percent  to  the  radiated  seismic 
energy. 

With  the  descent  of  the  isotherms  inside  the  slab,  could 
the  material  fail  in  a  mode  seismically  equivalent  Lc 
brittle  fracture  at  a  depth  of  600  km  or  more?  In  the  absence 
of  direct  laboratory  data  on  rock  failure  at  200  kb  confining 
pressure,  this  problem  can  only  be  explored  indirectly.  If 
we  assume  diffusion  creep  to  be  the  mechanism  of  deformation 
(Gordon  1965;  Orowan  1967;  Carter  &  Ave 'Lallemant  1970; 

Raleigh  &  Kirby  1970;  Weertman  1970;  Goetze  1971) ,  we  can 
calculate  the  strain  rate  e: 


e  =  5.37 


x  10*®  x  exp 


H+V  P 
ac 


RT 


.2.3 


(10) 


where  H  =  activation  enthalpy,  V  =  activation  volume,  P  = 

dC 

hydrostatic  pressure,  R  =  gas  constant,  T  =  temperature  (°K) , 
and  o  =  maximum  shear  stress.  We  can  also  compute  the 
effective  viscosity  n: 


e 


Taking  our  calculated  value  for  stress  (i.e.,  o  =  500  bars), 
and  H  =  100  kcal'mole,  V  =  40  cm3/mole  (Goetze  1971,  and 
personal  communication)  and  *ho  temperature  in  the  center 
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of  the  slab  from  Fig.  3,  we  find  within  20  km  of  the  center 

27 

the  effective  viscosity  is  greater  than  n  =  10  poise  for 
virtually  the  full  length  of  the  slab.  This  value  is  indeed 
very  high,  and  it  implies  a  very  low  creep  rate  inside  the 
slab.  The  slab  behaves  nearly  elastically,  and  uniform 
creep  cannot  relieve  the  stresses.  It  is  necessary  to  have 
either  accelerated  creep  or  brittle  fracture  to  occur  in 
response  to  stress  concentrations.  In  the  case  of  accel¬ 
erated  creep,  however,  the  maximum  deformation  must  take 
place  over  a  very  narrow  (low  aspect  ratio)  volume  to  produce 
the  equivalent  of  a  dislocation  model  for  the  seismic  source. 
This  mechanism  is  in  agreement  with  the  creep  model  described 
by  Griggs  (1972)  based  on  experiments  with  dunite. 

The  maximum  depth  of  earthquakes  is  also  of  great  interest. 
No  earthquakes  have  been  observed  below  a  depth  of  730  km, 
and  many  seismic  zones  terminate  at  shallower  depths.  The 
absence  of  wel 1-def i ned  travel  time  anomalies  for  deep  focus 
earthquakes  indicates  that  a  well-defined  slab  does  not  exist 
for  a  significant  distance  below  the  earthquakes  (Toksttz  et  al. 
1971;  Mitronovas  (<  I  sacks  1971;  Sen  Gupta  &  Julian  1973).  A 
broad  region  of  somewhat  reduced  temperature  or  a  weak  con¬ 
tinuation  of  the  slab  below  the  depth  of  maximum  earthquakes 
cannot  be  excluded  by  the  data. 

This  limited  depth  of  earthquakes  may  be  explained  by 
•  'jthei  the  absence  of  stresses  nece^Sc.ry  to  produce  earthquakes 
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or  by  the  absence  of  material  which  would  undergo  brittle 
failure.  These  two  are  related  since  the  temperature  plays 
a  strong  role  in  both  cases.  Below  a  depth  of  650  km, 
the  interior  of  the  descending  slab  rapidly  reaches  thermal 
equilibrium  because  of  efficient  radiative  heat  transfer. 

Yet  equilibration  through  thermal  conduction  is  not  likely 
to  be  the  sole  cause  of  a  maximum  depth  for  earthquakes,  as 
earthquakes  in  some  seismic  zones  such  as  Tonga  do  not  become 
progressively  rarer  with  depth  (Isacks  et  al.  1968) .  It  is 
likely  that  the  650  km  phase  change  rises  to  shallower  depths 
in  the  slab  and  produces  the  stress  for  earthquakes  below 
650  km  (Smith  &  Toksfiz  197?) .  The  seismicity  of  the  small 
deteched  segments  of  the  slabs  under  New  Hebrides  and  under 
Tonga  can  also  be  explained  by  the  stresses  due  to  this 
phase  boundary  and  the  resistance  of  the  lower  mantle. 

Mechanisms  of  Shallow  Earthquakes  in  Convergence  Zones 

There  are  a  large  number  of  shallow  earthquakes  that 
also  occur  in  the  convergence  zones,  under  or  near  the  trenches 
There  is  a  general  agreement  regarding  the  mechanisms  of 
these  shallow  earthquakes.  On  the  basis  of  locations  and 
source  mechanisms  these  can  be  divided  into  two  categories: 

1.  The  shallow  earthquakes  associated  with  normal 
faults  in  the  oceanic  side  of  the  trenches  are  probably 
due  to  tension  in  the  outer  layer  of  the  lithosphere  arising 
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from  the  bending  of  the  slab  and  the  upward  buckling 
(Stauder  1968b;  Hanks  1971;  Kanamori  1971)  .  Generally 
these  earthquakes  are  not  very  large  and  do  not  cut  across 
the  Lithosphere  (the  Sanriku  earthquake  of  1933  may  be  an 
exception) .  The  net  stress  in  the  lithosphere  may  be 
compressive  as  evidenced  from  the  earthquakes  discussed  below. 

2.  The  second  group  of  shallow  earthquakes  associated 
with  subduction  occurs  along  the  shear  zone  between  the  down¬ 
going  slab  and  the  continental  lithosphere.  Many  major 
earthquakes  are  associated  with  thrust  zones  and  fall  in  this 
category.  The  dip  of  the  fault  plane  is  generally  less  than 
about  30° .  Some  examples  of  these  are  the  Chile  earthquake 
of  1960,  the  Alaska  earthquake  of  1964  (Plafker  1972),  and 
the  Rat  Island  earthquake  of  1965  (Stauder  1968a) .  This 
type  of  earthquake  is  also  prevalent  in  regions  where  con¬ 
tinental  lithospheric  plates  converge.  The  earthquakes  in 
the  southern  Iran  thrust  belt  (Nowroozi  1972)  are  good 
examples  of  this  category. 

We  should  clarify  that  in  all  the  above  areas  some  events 
not  matching  the  two  generalized  categories  do  occur  (e.g., 
strike  slip  faults).  Considering  the  immensely  complicated 
crustal  features,  tears  in  the  plates,  and  localized  stress 
concentrations,  faults  along  other  orientations  and  earth¬ 
quakes  of  differing  source  mechanisms  would  and  do  occur. 
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Thus  the  previous  generalizations  about  earthquake  charac¬ 
teristics  represent  the  predominant  types  rather  than  all 
earthquakes  in  the  convergence  areas. 
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CONCLUSIONS 

Our  understanding  of  deep  earthquakes  and  the  evolution 
of  the  descending  lithosphere  relies  upon  what  inferences 
can  be  extracted  from  observations  and  from  predictions  using 
theoretical  and  laboratory  work.  Here  theoretical  calcu¬ 
lations  of  the  thermal  regime,  the  seismic  velocity  structure, 
and  the  stress  in  the  descending  lithosphere  provide  the 
predictions;  and  these  are  in  good  agreement  with  relevant 
observations.  A  major  contribution  of  the  theoretical 
calculations  -  in  addition  to  verifying  the  hypothetical 
models  -  has  been  their  ability  to  identify  the  significant 
parameters  that  control  the  evolution  and  properties  of  the 
descending  lithosphere.  Some  major  conclusions  of  this 
study  are: 

1.  The  temperatures  inside  the  slab  are  strongly 
controlled  by  the  time  elapsed  after  initiation  of 
the  descent  and  by  the  conductivity,  especially  the 
radiative  term.  The  gross  structure  of  the  thermal 
regime  is  insensitive  to  reasonable  variations  in  the 
parameters . 

2.  The  maximum  depth  of  penetration  is  strongly  influenced 
by  the  650  km  phase  change,  and  by  the  increase  in  the 
radiative  contribution  to  effective  thermal  conductivity. 
Temperature  anomalies  do  not  extend  beyond  700  km  if  the 
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above  phase  change  is  exothermic.  In  all  cases, 
however,  the  slab  is  thermally  assimilated  into  the 
mantle  at  some  depth  shallower  than  about  800  km. 

3.  For  convergence  of  continental  plates,  the  underthrusting 
^  thosphere  generally  does  not  penetrate  very  deep. 
Orogeny  associated  with  continental  collisions  most 
likely  occurs  either  during  the  collision  or  later 

when  the  subducted  continent  becomes  heated. 

4.  From  the  comparison  of  travel  time  anomalies  associated 
with  the  LONGSHOT  explosion  and  the  calculated  values 
based  on  theoretical  temperatures  and  velocities,  the 
intermediate  depth  earthquakes  are  located  within  the 
coolest  region  of  the  slab  under  the  Aleutians.  This 
also  seems  to  be  the  case  in  the  Tonga  region  and  under 
Japan  as  inferred  from  the  velocity  anomalies  and 
attenuation  of  seismic  waves. 

5.  Stress  calculations  shov/  that  the  maximum  shear  stress 
occurs  within  the  coolest  region  of  the  slab  for  a 
reasonable  temperature-dependent  rheology. 

6.  The  regional  stress  patterns  depend  upon  the  temperature 
field,  the  slab  rheology,  dip  of  the  slab,  and  mantle 
support.  The  last  item,  the  mantle's  support  of  the 
slab,  may  be  influenced  by  the  descent  velocity, 
convection,  phase  changes,  and  lateral  differences 
within  the  mantle  such  as  under  oceans  or  continents. 
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7.  The  magnitude  of  the  calculated  shear  stress  within 

the  slab  is  about  500-1000  bars.  This  is  a  reasonable 
value  as  inferred  from  estimates  of  the  stress  drop 
for  deep-focus  earthquakes.  The  orier nation  of  the 
theoretical  principal  stresses  are  in  very  good 
agreement  with  those  determined  from  focal  mechanism 
studies  of:  intermediate  and  deep-focus  earthquakes. 

The  center  of  the  slab  is  sufficiently  cold  to  allow 
either  accelerated  creep  or  brittle  fracture  to  occur 
in  response  to  stress  concentrations. 
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FIGURE  CAPTIONS 

1  Schematic  representation  of  the  computational 

scheme,  input  parameters,  and  boundary  conditions 
used  to  calculate  thermal  models  of  slabs  (Modi¬ 
fied  after  Minear  &  Toksttz  1970a) . 

Fig.  2  Temperatures  inside  the  base  slab  model  after  3.6 

m.y.  and  7.1  m.y.  were  calculated  using  the  parameters 
given  in  the  text.  The  velocity  of  the  slab  is 
8  cm/yr . 

Fig.  3  The  slab  in  Fig.  2  after  10.7  m.y.  Penetration 

through  the  phase  transition  at  600  km  depth,  and 
the  latent  heat  of  the  pnase  transition  help  heat 
and  thermally  equilibrate  the  slab  interior  below 
this  depth. 

Fig.  4  Temperature  field  of  the  slab  calculated  using 

Schatz  i  Simmons  (1972)  conductivity.  All  other 
parameters  same  as  those  in  Fig.  3.  Due  to  lower 
conductivity,  this  slab  penetrates  the  600  km 
phase  change. 

Fig.  5  Temperature  field  of  the  slab  (shown  in  Fig.  3) 
after  16  m.y.  The  penetration  of  the  slab  was 
stopped  at  10.7  m.y.,  and  the  slab  remained  at  rest 
fot  5.3  m.y. 
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Fig  .  6 


Fig .  7 


Fig .  8 


Fig .  9 


The  thermal  fields  of  the  slab  in  Fig.  5,  after 
additional  rest  periods.  The  temperature  anomaly 
becomes  progressively  broader  with  time. 

Temperature  field  in  a  slab  penetrating  into  a 
mantle  having  a  "convective"  geotherm.  The  param¬ 
eters  used  for  this  slab  are  self-consistent  with 
the  geotherm.  Penetration  of  the  slab  was  stopped 
at  10.7  m.y.,  and  the  lower  right  siab  was  allowed 
to  equilibrate  with  the  mantle  for  3.6  m.y.  The 
phase  boundary  at  650  km  was  not  taken  into  account 
in  tnis  model,  but  the  higher  phase  (350  km)  boundary 
was  included. 

Temperature  field  in  subducted  continent  after 
10  m.y.  with  extreme  model  of  frictional  heating 
along  the  fault  plane.  The  highly  radioactive 
continental  crust  is  indicated  by  hatching. 

The  temperature  in  a  subducted  continent  after 
7.7  m.y.  (above),  and  after  remaining  stationary 
for  an  additional  47  m.y.  (below).  A  thermal 
maximum  has  formed  in  the  subducted  crust.  "M" 
indicates  the  base  of  the  crust. 


i 


-191- 


Fig.  10  Observed  and  calculated  P-wave  travel  time  anomalies 
for  LONGSHOT  explosion.  The  observed  delays  are 
station  corrected  and  averaged  over  intervals  of 
azimuth  and  distance  as  given  by  Abe  (1972a) .  The 
fit  is  obtained  by  assuming  the  strike  of  the  arc 
is  N82W,  and  adjusting  the  slab  relative  to  the 
LONGSHOT  location. 

Fig.  11  Temperature  regime  and  the  distribution  of  earth¬ 
quake  hypocenters  in  the  Central  Aleutians  near 
Amchitka.  The  thermal  model  is  obtained  by  combining 
the  isotherms  from  two  slabs  dipping  at  different 
angles.  The  earthquake  hypocenters  are  from 
Engdahl  (1971) .  The  focal  mechanisms  of  two  repre¬ 
sentative  earthquakes  show  normal  faulting  in  the 
Oceanside  of  the  trench  and  thrusting  in  the  con¬ 
vergence  zone  beneath  the  arc.  Note  how  intermediate 
earthquakes  are  in  the  coolest  region  of  the  slab 
and  not  in  the  shear  zone. 

Fig.  12  Density  anomaly  in  gm/cm'*  within  tne  descending  slab 


using  a  typical  thermal  regime  (in  this  case  Fig.  2 
of  Smith  &  Toksflz ,  1972).  Only  thermal  volume 
contraction  .lb  used  for  the  density  calculations. 
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Fig.  13  Calculated  maximum  shear  stress  models  for  two 

different  mantle  supports,  using  the  30°  dipping 
slab  shown  in  Fig.  12.  On  the  left  the  mantle 
support  is  given  in  terms  of  an  'effective'  Young's 
modulus,  E,  versus  depth.  Section  A-A'  depicts  the 
modulus  within  the  slab  which  is  temperature 
dependent.  The  arrows  indicate  the  direction  of  the 
principal  stress  closest  to  down-dip,  and  the 
contours  are  maximum  shear  stress  in  bars.  Shear 
stresses  within  the  mantle  are  not  contoured,  but 
are  generally  lower  than  200  bars. 

Fig.  14  Comparison  of  theoretical  model  (top  model.  Fig.  13) 
to  the  seismicity  of  Japan  given  by  Utsu  (1971) . 

Fig.  15  Theoretical  stresses  and  earthquake  focal  mechanisms 
for  Honshu  and  Kermadec .  On  the  left  are  the 
theoretical  models.  Both  models  have  increasing 
support  with  depth  for  the  mantle  and  strong 
temperature  dependence  for  the  slab  rheology,  although 
each  is  slightly  different.  The  dips  of  the  slabs 
are  specified  by  the  seismicity.  On  the  right  are 
focal  mechanism  solutions:  for  each  region,  taken 
from  Isacks  &  Molnar  (1971)  ,  superimposed  upon 
the  plane  of  seismicity.  Open  circles  show  down-dip 
compression;  filled  circles  indicate  down-dip  extension. 
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4 .  ARRAY  STUDIES 

4 • 1  Automatic  Classification  of  Seismic  Detections 

From  Large  Aperture  Seismic  Arrays  by  S.  Shlien 
(Abstract) 

The  large-aperture  seismic  arrays  in  Montana  (LASA) 
and  Norway  (NORSAR)  make  on-line  signal  processing  a 
necessity  if  these  arrays  are  to  be  used  at  their  full 
capability.  Using  the  outputs  of  the  detection  processors 
of  the  respective  arrays,  the  feasibility  of  automatic 
classification  of  seismic  signals  into  the  various  body 
phases  P,  PKP ,  PcP ,  ScP ,  SKP,  PP ,  PKKP  and  P'P'  was 
confirmed.  It  was  shown  how  these  later  phases  can  be 
used  to  advantage  in  improving  the  location  capability 
using  the  combination  of  the  two  arrays. 

One  of  the  byproducts  of  this  study  was  an  estimation 
of  the  detection  and  location  capabilities  of  the  arrays. 
It  was  estimated  that  LASA  detects  more  than  50  real 
seismic  signals  a  day,  of  which  less  than  10%  are  due  to 
later  phases.  LASA ' s  detection  capability  extends  almost 
one  body  wave  magnitude  below  ERL's  capability  based  on 
reported  epicenters .  The  discrimination  between  very 
weak  seismic  signals  and  false  alarms  due  to  spurious 
noise  was  found  difficult  on  the  basis  of  only  the 


-209- 


detection  logs. 

Only  a  little  more  than  8  earthquakes  a  day  were 
found  common  between  LASA  and  NORSAR  arrays.  It  is 
expected  that  this  number  will  increase  with  the  im¬ 
proved  signal  processing  that  the  two  arrays  recently 
implemented. 

4 . 2  Automatic  Event  Detection  and  Location  Capabilities 
of  Large  Aperture  Seismic  Arrays  by  S.  Shlien  and 
M.N.  ToksBz 
Summary 

The  detection  and  location  capabilities  of  large 
aperture  seismic  arrays  (LASA  in  Montana  and  NORSAR  in 
Norway)  were  determined  on  the  basis  of  the  Detection 
Logs,  Summary  Bulletins  and  NOAA-ERL  Preliminary  Determin¬ 
ation  of  Epicenters.  Estimates  of  signal  and  false  alarm 
distributions  indicated  tnat  LASA  can  detect  on  the 
average  60  seismic  signals  a  day  and  NORSAR  about  20.  The 
probability  of  detection  of  an  event  was  estimated  for 
LASA  and  NORSAR  as  a  function  of  distance  and  magnitude. 

A  method  of  obtaining  theoretical  estimates  of  the 
location  capability  of  an  array  as  a  function  of  beam 
resolution,  phase  type,  and  slowness,  dT/dA,  is  described. 
Theoretical  estimates  of  the  array's  location  capability 
compared  to  observed  estimates  imply  that  the  arrays  are 
operating  close  to  their  expected  limits. 
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INTRODUCTION 

The  detection  and  location  of  earthquakes  in  real  time 
can  be  accomplished  with  the  utilization  of  large  aperture 
seismic  arrays .  About  30  earthquakes  a  day  are  reported 
in  the  LASA  Summary  Bulletin  on  a  routine  basis.  In  this 
paper  the  present  detection  and  location  capabilities  of 
LASA  and  NORSAR  were  estimated.  Because  of  improvements  of 
processing  schemes,  the  present  paper  updates  some  earlier 
studies  (Dean  et  al . ,  1971;  Evernden,  1971;  IBM  Final  Report, 
1972) . 

The  preparation  of  the  Summary  Bulletins  can  be  divided 
into  three  basic  steps,  indicated  by  the  b^ock  diagram  in 
Figure  1.  In  the  first  step,  the  Detection  Processor  attempts 
to  flag  all  the  possible  seismic  signals  arriving  at  the 
array  of  short-period  seismometers.  These  potential  signals 
are  recorded  into  the  Detection  Log  in  real  time  together 
with  rough  estimates  of  the  Signal  to  Noise  Ratio  (SNR) , 
signal  slowness  vector,  and  signal  intensity.  In  the  following 
step,  the  Event  Processor  selects  the  detections  with  the 
larger  SNR  and  then  refines  the  estimates  of  signal  slowness 
and  intensity  of  these  detections.  The  event  epicenter, 
origin,  and  magnitude  are  computed  from  the  estimates. 

The  last  step  involves  the  editing  of  the  output  of  the 
Event  Processor  into  the  Summary  Bulletin.  This  step 
is  performed  by  analysts  who  verify  the  output  of 


the  Event  Processor  using  the  actual  waveforms  from  the 
subarrays . 


In  this  study  both  the  Detection  Logs  and  Summary 
Bulletins  from  LASA  and  NORSAR  were  used  to  evaluate  the 
capabilities  of  the  arrays.  Data  were  available  for  two 
periods:  the  "summer  period"  from  Kay  to  mid-August  1971  and 
the  "winter  period"  of  20  February  to  19  March  1972.  The 
epicenter  determinations  published  by  the  Environmental 
R«ea£Ch_Laboratori  (ERL)  were  used  as  an  outside  standard. 

Since  ERL  uses  the  LASA  bulletin  as  a  primary  source  in  addition 
to  other  single  seismic  stations,  ERL  estimates  have  a  small 
dependence  on  LASA  measurements.  Except  for  magnitude 
estimates  of  small  earthquakes  this  effect  may  be  neglected. 

Heavy  emphasis  was  placed  upon  the  Detection  Logs  in  the 
estimation  of  the  detection  capabilities  of  the  arrays.  The 
Detection  Logs  were  the  most  complete  list  of  possible  signals 
observed  at  the  arrays.  A  detailed  description  of  these  logs 
is  given  in  the  second  section.  The  location  capabilities  of 
the  arrays  were  determined  mainly  from  the  Summary  Bulletins  which 
contained  the  most  accurate  estimates  of  the  signal  paras, eters. 
Those  are  discussed  in  the  fourth  section. 


THE  DETECTION  PROCESSOR  AND  THE  DETECTION  LOG 
The  Detection  Processor  for  a  large  aperture  seismic 
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array  was  designed  and  developed  by  IBM  and  went  into  oper¬ 
ation  as  of  April  1969.  Except  for  a  few  minor  details  the 
LASA  and  NORSAR  Detection  Processors  run  in  a  similar  manner. 
Descriptions  of  the  LASA  and  NORSAR  systems  are  given  in  the 
IBM  Final  Report  (1972)  and  in  Bungum  et  al^.  (1971)  ,  respec¬ 
tively.  Here  we  include  a  brief  summary  of  the  procedures, 
in  order  to  facilitate  the  understanding  of  the  later  parts 
of  the  paper. 

The  Detection  Processor  accomplishes  two  basic  oper¬ 
ations:  (1)  it  increases  the  SNR  by  beam  forming  and  fil¬ 

tering,  and  (2)  it  decides  whether  a  signal  is  present.  A 
teieseismie  signal  arrives  at  the  array  approximately  in  the 
form  of  a  plane  wave  with  a  specific  slowness  and  azimuth, 
depending  upon  the  location  of  the  earthquake  and  the  phase 
type.  The  noise,  on  the  other  hand,  is  generally  diffuse 
and  incoherent.  Beam  forming  increases  the  SNR  by  screening 
out  the  noise  from  all  directions  and  slowness  except  that  of 
the  actual  signal,  while  leaving  the  signal  approximately 
intact.  Several  hundreds  of  beams  are  generated  in  real  time 
by  the  LASA  and  NORSAR  Detection  Processors  by  delaying  and 
summing  the  subarray  signals.  Each  beam  is  presteered  towards 
a  specific  direction. 

The  Detection  Processors  deploy  two  sets  of  overlapping 
beams.  The  first  partition  consists  of  narrow  beams  which 
cover  mainly  the  seismic  regions  and  other  areas  of  special 
interest.  This  partition  is  plotted  on  the  wor Ld  map  in 
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Figure  2  for  the  P  and  PKP  phases.  The  second  partition 
consists  of  broader  beams  which  cover  the  entire  signal 
space . 

Each  of  the  beams  is  filtered  and  rectified  by  the  De¬ 
tection  Processor.  The  filter  was  designed  to  deemphasize 
those  frequency  components  where  the  SNR  is  low.  In  the 
case  of  the  LASA  array  the  signal  is  confined  tc  a  narrow 
band  centered  at  1  Hz.  The  signal  at  NORSAR  covers  a  broader 
frequency  band,  between  .5  and  5.  Hz.  The  rectified  beams 
then  pass  through  two  integrators  of  different  time  durations 
These  integrators  compute  a  Short  Term  Average  (STA)  and  a 
Long  Term  Average  (LTA) .  The  LTA  is  determined  over  a  28- 
second  interval  and  is  supposed  to  be  a  measure  of  the 
natural  noise.  The  STA  is  computed  for  a  1.8-second  time 
interval  and  is  a  measure  of  the  amount  of  signal,  if  present 
The  LTA  and  STA  are  measured  in  so-called  quantum  units  where 
1  quantum  unit  was  set  at  .0028  millimicrons  for  LASA  and  at 
.0010  millimicrons  for  NORSAR.  The  STA  and  LTA  are  updated 
every  0.6  and  1.8  seconds  respectively.  If  20  log^g 
(STA/LTA)  is  above  8  db  for  at  least  two  seconds,  then  the 
particular  beam  is  declared  to  be  in  the  detection  state. 

A  large  signal  will  usually  trigger  several  beams 
simultaneously.  The  beams  with  the  maximum  STA  in  each  of 
the  beam  partitions  are  recorded  onto  the  detection  log  for 
that  particular  time  cycle.  A  large  seismic  signal  usually 
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has  several  bursts  of  energy  so  that  as  many  as  15  beam 
detections  could  be  recorded  for  just  a  P  phase. 

The  detection  log  is  a  record  of  the  parameters  of 
every  detection  declared  by  the  processor.  The  beam  number, 
the  Maximum  STA  (MSTA) ,  the  LTA  prior  to  the  detection  state 
and  the  start  and  stop  times  of  the  detection  are  recorded 
into  the  log  in  chronological  order.  On  the  average  there 
were  450  detection  entries  in  the  LASA  log  and  140  detec¬ 
tion  entries  in  the  NORSAR  log  per  day  for  the  period  20 
February  to  19  March  1972.  Most  of  these  detections  occured 
in  groups  of  two  or  more.  The  cumulative  number  of  detec¬ 
tion  groups  per  day  versus  the  maximum  MSTA  of  the  group 
is  plotted  for  LASA  and  NORSAR  for  the  winter  period  in  the 
upper  half  of  Fiy'.re  3.  (A  detection  was  defined  to  belong 
to  the  same  group  if  it  occurred  within  30  seconds  of  the 
previous  detection) . 

The  detection  groups  with  large  MSTA  were  confined 
mostly  to  beams  pointed  towards  seismically  active  areas. 

In  contrast  the  numerous  detection  groups  with  low  MSTA's 
occurred  indiscriminately  among  the  beams  directed  towards 
seismic  and  aseismic  areas.  This  observation  coupled  with 
the  fact  that  detections  with  low  MSTA  have  a  low  SNR, 
strongly  suggested  that  the  low  MSTA  detections  are  pre¬ 
dominantly  false  alarms. 
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Estimation  of  the  distribution  of  false  alarms  and 
signals  could  be  determined  only  indirectly.  The  top  part 
of  Figure  3  shows  the  cumulative  distribution  of  detec¬ 
tions  versus  the  maximum  MSTA  of  the  group.  For  large 
MSTA's  the  SNR  precludes  the  possibility  of  the  detection 
group  being  a  false  alarm  and  the  MSTA  distribution  obeys 
a  Pareto  (power  law)  distribution.  This  distribution  was 
observed  by  Milne  and  Davenport  (1969)  and  follows  directly 
from  Gutenberg  and  Richter's  magnitude  distributions.  On 
account  of  the  clipping  of  very  large  signals  by  the  Detection 
Processor,  the  distribution  departs  from  the  power  law  at 
high  MSTA.  As  the  MSTA  decreases  below  a  certain  value  the 
number  of  detections  starts  to  increase  abnormally  due  to  the 
inclusion  of  many  false  alarm  groups.  The  MSTA  distribution 
of  the  signals  was  estimated  for  low  MSTA  by  extrapolating 
the  distribution  to  the  right  of  MSTA  =  300  for  LASA  and 
MSTA  —  2000  for  NORSAR  with  a  function  of  the  form 

N  =  A (MSTA) “b  (1) 

where  A  and  b  were  chosen  to  yield  the  best  fit.  This  is 
indicated  by  the  dashed  line. 

The  false  alarm  distribution  was  determined  by  subtracting 
the  presumed  cumulative  signal  distribution  from  the 
cumulative  detection  group  distribution.  The  probability 
of  a  detection  being  a  false  alarm,  given  MSTA,  was  determined 
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by  taking  the  ratio  of  false  alarm  groups  to  detection 
groups  in  given  MSTA  ranges.  This  was  plotted  in  the  lower 
half  of  Figure  3.  The  false  alarm  probability  functions 
obtained  in  this  manner  agreed  with  those  determined  inde¬ 
pendently  using  a  set  of  aseismic  beams  described  by  Shlien 
(1972)  . 

It  is  apparent  that  in  order  to  limit  false  alarm 
occurrences  the  MSTA  threshold  must  be  raised  above 

200  for  LASA  and  750  for  NORSAR.  With  MSTA  =  200,  as  many 
as  50  seismic  signals  are  detected  in  the  case  of  LASA. 

In  the  next  section  the  detection  capabilities  of 
LASA  and  NORSAR  are  estimated  as  a  function  of  magnitude  and 
distance  of  the  event. 


DETECTION  CAPABILITIES  OF  LASA  AND  NORSAR 

The  detection  capabilities  of  seismic  instruments  are 
limited  primarily  by  microsei smic  noise.  The  noise  level 
is  considerably  larger  at  NORSAR  in  the  signal  band  on 
account  of  its  proximity  to  the  coast  (IBM  Final  Technical 
Report,  1971;  Felix  e_t  al_.  ,  1972)  .  A  large  aperture  seismic 
3^^«*y  allows  one  to  enhance  the  SNR  by  appropriate  beam 
forming.  If  the  noise  is  uncorrclated  and  the  signal  is 

S' 

coherent  between  the  sensors,  the  maximum  St\'R  gain  achievable 

wj  th  an  array  joes 
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as  /~N  where  N  is  the  total  number  of  sensors. 

For  two  reasons  neither  LASA  nor  NORSAR  can  attain 
this  optimum  gain.  First,  the  signal  amplitude  and  shape 
vary  by  almost  an  order  of  mignitude  across  the  array.  This 
is  believed  to  be  caused  by  path  differences  between  source 
and  sensors  and  by  the  multipathiny  and  scattering  of  the 
signal  underneath  the  array  (Mack,  1969;  Lamer,  1970)  . 

The  signal  variations  are  consistent  with  events  from  the 
same  location  but  change  unpredic tably  when  the  epicenter 
is  shifted  by  a  few  degrees.  Finally,  the  seismic  signal 
arrives  as  a  distorted  plane  wave  due  to  local  and  deep 
lateral  heterogeneities  (Greenfield  and  Sheppard,  1969? 

Davies  and  Sheppard,  1972) .  Unless  special  velocity-dependent 
station  corrections  were  incorporated  into  the  beam  forming 
process  severe  signal  losses  would  be  sustained.  Both  LASA 
and  NORSAR  had  an  improved  set  of  station  corrections  since 
February  1972.  As  a  result  there  has  been  a  slight  improvement 
in  the  detection  and  location  capabilities  of  LASA  and  NORSAR 
since  that  time. 

Using  the  ERL  preliminary  epicenter  determinations 
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as  a  standard  we  estimated  the  capabilities  of  LASA  and 
NORSAR  by  counting  the  number  of  matches  that  could  be 
made  with  the  detection  log.  The  arrival  time  of  the  P 
or  PKP  phase  was  calculated  at  either  array  using  the 
ERL  catalog  and  travel  time  tables.  The  signal  was 
taken  to  be  detected  by  the  array  if  it  could  be  matched 
to  a  detection  in  the  log  within  ±10  seconds  of  the  ex¬ 
pected  arrival  time.  Periods  when  the  Detection  Processor 
was  down  were  taken  into  account  in  the  analysis. 

To  reduce  the  probability  of  a  spurious  match  with  a  signal  or 
false  alarm,  detections  below  MSTA  =  250  for  LASA  and 
MSTA  =  1000  for  NORSAR  were  disregarded.  With  this  rule 
it  was  estimated  that  the  probability  of  making  a  false 
match  was  about  1  percent  for  either  array.  In 
Figure  4  the  percentage  of  ERL  events  detected  per  10 
degree  distance  intervals  are  plotted  for  LASA  and 
NORSAR.  Data  for  the  20  February  to  19  March  1972  period 
are  indicated  by  the  dashed  lines.  Both  arrays  had  made 
slight  improvement  in  their  detection  capabilities. 

Most  of  the  improvement  is  due  to  the  incorporation  of 
more  accurate  station  corrections  into  the  systems  be¬ 
fore  the  winter  period  data  were  taken. 

The  detection  capabilities  of  the  arrays  are  con¬ 
sistently  inferior  for  local  events  and  dist.ant  events 
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in  the  vicinity  of  the  shadow  zone  (100-110  degrees). 

Local  events  are  difficult  to  detect  by  the  present  pro¬ 
cessors  for  the  following  three  reasons.  (1)  The  signal 
is  emergent  and  spread  out,  (2)  the  wavefront  is  too  spher¬ 
ical  and  cannot  be  approximated  by  a  plane  wave,  and  (3) 

the  dT/dA  varies  so  rapidly  with  distance  in  this  range 

LASA  has  no  beams  in  the  near  zone, 
that  a  prohibitive  number  of  beams  would  be  needed.  Alie”ond 

90  degrees  the  P  wave  becomes  diffracted  and  highly  atten¬ 
uated  by  the  core-mantle  boundary.  Low  magnitude  earth¬ 
quakes  are  easily  missed.  In  the  diffraction  zone  the 
dT/dA  versus  A  curve  becomes  nearly  constant.  Hence  it 
is  very  difficult  to  locate  earthquakes  from  this  zone. 

The  LASA  Summary  Bulletin  does  not  attempt  to  report  any 
events  beyond  100  degrees.  NORSAR  appears  to  miss  an  un¬ 
usual  number  of  events  coming  from  the  Aleutians  and  the 
western  part  of  the  United  States.  This  is  evidenced  by 
the  dip  in  the  detection  capability  beyond  60  degrees. 

In  Figure  5  the  detection  capabilities  of  LASA  and 
NORSAR  are  shown  as  a  function  of  magnitude.  The  percents 
of  ERL  events  detected  between  0  and  80  degrees  (upper  half) 
and  between  80  and  180  degrees  (lower  half)  distance  from 
the  array  in  question  were  plotted  versus  magnitude  units. 

Newer  data  from  the  winter  period  (1972)  are  distinguished 
by  the  dashed  lines.  Due  to  the  rapidly  decreasing  sample 
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sizes  at  the  higher  magnitudes,  the  percent  of  events  de¬ 
tected  sometimes  decreases  with  increasing  magnitude.  LASA 
detects  almost  all  ERL  events  greater  than  magnitude  4.5 
and  less  than  80  degrees.  Considerable  improvement  in  both 
NORSAR's  and  LASA's  detection  capability  after  February  1972 
is  evident  in  the  lower  plots. 

It  is  evident  that  NORSAR's  detection  capability  is 
overall  inferior  to  LASA's.  Unlike  LASA,  a  considerable 
amount  of  the  signal  energy  in  NORSAR  is  concentrated  at 
high  frequencies  (Harley,  1972).  However,  since  the  noise  is 
also  higher  in  the  high  frequency  band  at  NORSAR  the  problem 
of  discriminating  the  signal  from  noise  is  more  acute.  In 
addition,  the  effect  of  inaccurate  station  corrections  is 
worse  for  high  frequency  signals. 

The  analysis  so  far  has  shown  that  LASA  can  detect 
over  90  percent  of  the  reported  (ERL)  events  in  the  distance 
range  of  30  to  80  degrees  and  that  NORSAR  likewise  can 
detect  close  to  60  percent  of  the  ERL  events.  It  is  not 
apparent  so  far  how  many  additional  events  the  arrays  de¬ 
tect  that  are  not  on  the  ERL  list.  To  make  this  determination 
we  resorted  to  the  Summary  Bulletins  generated  by  LASA  and 
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In  Figure  6  the  frequency-magnitude  distribution  of 
events  reported  in  the  Summary  Bulletins  are  plotted  with 
the  ERL  frequency-magnitude  distribution.  Events  further 
than  95  degrees  of  the  array  in  question  were  not  included 
in  either  distribution.  It  is  seen  that  LASA  reports  many 
more  low  magnitude  events  than  ERL.  The  peak  of  the  LASA 
distribution  m^  =  3.7  is  one  magnitude  unit  to  the  left  of 
the  ERL  distribution.  NORSAR  apparently  reports  slightly 
more  low  magnitude  events  than  ERL  and  the  peak  lies  0.5 
magnitude  units  to  the  left  of  ERL's  mode. 

The  leftward  shift  of  the  frequency-magnitude  distri¬ 
bution  can  also  be  partially  attributed  to  a  bias  in  the 
magnitude  estimates.  Correlation  plots  of  LASA  and  NORSAR' s 
magnitudes  versus  ERL  magnitude  are  shown  in  Figure  7.  They 
suggest  that  part  of  the  shift  in  the  distributions  can  be 
explained  by  small  magnitude-dependent  biases. 

These  results  have  demonstrated  that  LASA  and  NORSAR 
have  detection  capabilities  comparable  to  ERL  at  least  for 
earthquakes  in  the  range  of  30  to  80  degrees  distance. 

Similar  results  have  been  disclosed  by  the  IBM  Final  Technical 
Report  (1971)  and  by  Dean  et  al.  (1971) .  On  the  other  hand 
the  location  capability  of  an  array  would  not  be  as  accurate 
as  that  obtained  with  a  network  of  seismometers  due  to  its 
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finite  aperture.  Epicenter  locations  accurate  to  a  tenth 
of  a  degree  can  be  achieved  with  seismic  stations  suitably 
spaced  around  the  earthquake's  focus  as  was  shown  by 
Evernden  (1971) . 

In  the  next  section  a  theoretical  analysis  of  the 
location  capability  of  the  arrays  was  carried  out  and  the 
actual  performance  of  arrays  in  locating  events  is 
described. 

LOCATION  CAPABILITY  OF  LASA  AND  NORSAR 

With  a  large  aperture  seismic  array  one  can  determine 
the  direction  and  velocity  of  the  arriving  signal.  This 
permits  an  estimate  of  the  event's  epicenter.  Studies  on 
the  location  capabilities  of  seismic  networks  and  arrays 
have  been  performed  by  Evernden  (1971)  and  Gjoystdal  et  al . , 
(1972)  . 

For  the  P  phase  there  is  a  single  valued  transformation 
between  the  signal's  inverse  phase  velocity  (slowness) 
vector  and  the  epicenter's  location.  The  distance.  A, 
of  the  event  depends  only  on  the  magnitude  of  the  slowness. 
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dT/d.  ,  and  the  azimuth  of  the  event  is  the  same  as  that 
of  the  signal. 

Due  to  the  limited  resolution  of  the  array  there  is 
always  an  uncertainty  a  in  the  determination  of  the  slow¬ 
ness.  The  size  of  this  uncertainty  is  proportional  to  T/A 
where  T  is  the  period  of  the  signal  and  A  is  the  aperture 
of  the  array.  The  associated  error  in  the  event's  location 
is  dependent  upon  its  distance  from  the  array  on  account  of 
the  nonlinear  transformations  involved. 

Distance  is  determined  from  dT/dA  using  an  empirical 
transformation  A  =  A (dT/dA)  based  on  the  present  models 
of  the  earth.  Since  each  phase  type  has  its  own  transfor¬ 
mation  it  is  important  to  have  the  signal  properly  classi¬ 
fied  at  this  point.  Supposing  the  measured  inverse  phase 
velocity  is  (dT/dA)  with  uncertainty  a,  then  the  un¬ 
certainty  in  the  corresponding  distance  is  given  by 


6A  =  A  [  (dT/dA)  ^  +  o/2]  -  A[(dT/dA)Q  -  o/2]  (2) 


If  o  is  sufficiently  small,  then  this  can  be  approximated 
by 


5A  = 


odA  (dT/dA) 
d (dT/dA)  I 

O' 


d2T 


dA2 


(3) 


o 


In  Figure  8,  dT/dA  and  d^T/dA^  curves  are  plotted 
versus  A  after  Herrin  et  al.  (1968)  for  the  P  phase. 
Around  90  degrees  d  T/dA  approaches  zero  rapidly,  re¬ 
sulting  in  very  poor  distance  determinations.  This  effect 
is  apparent  in  Figure  9, where  the  distance  and  azimuth  of 
all  earthquakes  triggering  specific  beams  in  the  LASA  high 
resolution  beam  partitions  over  a  certain  time  period  are 
plotted.  Though  the  beams  have  identical  resolution  in 
inverse  velocity  space,  it  is  very  clear  that  closer  to 
the  shadow  zone  the  region  of  epicenters  that  can  trigger 
the  same  beam  becomes  broader.  The  areas  represented  in 
the  figure  (Aleutians,  Japan,  and  Tonga)  are  extensive 
seismic  regions. 

The  azimuth  i esolving  power  of  an  array  is  also  de¬ 
pendent  upon  dT/c £ .  If  dT/dA  is  small  in  comparison  to 
o,  the  zone  of  unce'tainty  of  the  slowness  vector  is 
centered  close  to  the  origin  and  the  uncertainty  in  azimuth 
is  large.  This  corresponds  to  a  seismic  signal  coming 
nearly  vertical . 

In  order  to  es<  r.ate  analytically  the  error  in  azimuth 
it  was  assumed  that  the  error  in  the  slowness  vector  dT/dA 
determination  is  normally  distributed  with  zero  mean  and 
standard  deviation  o.  Let  u  be  the  magnitude  and  a 
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the  azimuth  of  the  actual  slowness  vector,  and  let  v  be 
the  measured  magnitude  and  3  the  measured  azimuth.  Then 
by  assumption  the  probability  of  measuring  v  and  3 
given  that  the  actual  values  are  u  and  a  is 


P  (v,  3  |  u,a)  =  ~^2  exp 


,22  2 

-[v  +u  -2vu  cos(a-3)]/2o 


(4) 


0  f  v  1  °° 

0  5  a  f  2u 


(This  is  the  same  model  for  the  radar  problem  of  narrow  band 
signal  with  additive  normal  noise.) 

Then 


p(3|u,a)  =  /p(v,3  |u,a)dv  (5) 

=  2rf  exP(_a0)  U  +  2/?aocosY1+,-rf^a.oCOSY)exp  (a^cos2y)  ] 
2  2  2 

where  aQ  =  u  /2o  and  y  =  8-a  (Schwartz,  1970).  For 
2 

large  aQ,  the  curve  may  be  approximated  by  a  Gaussian 
diStnoution  with  mean  a  and  variance  (2a^)-^.  The  stan¬ 
dard  error  of  8  is  plotted  versus  aQ  in  Figure  10.  (The 
standard  error  was  determined  numerically  for  small  a  .) 
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As  an  example,  the  theoretical  resolution  of  LASA  over 
a  frequency  band  1.0  -  2.0  Hz  is  about  a  =  0.2  seconds/degree 
(IBM  Final  Report,  1972) .  This  corresponds  to  an  error  of 
distance  of  6  A  =  2.5  degrees  in  the  distance  range  of  40 
to  85  degrees  and  an  error  of  6  A  =  20  degrees  in  the  shadow 
zone.  Clearly,  azimuth  error  is  much  smaller  than  distance 
error.  Azimuth  errors  do  not  become  appreciable  except  for 
phases  with  small  dT/dA  such  as  PKP  (DF  branch),  P'P'  (DF 
branch)  and  PKKP  (BC  branch) . 

This  analysis  had  not  considered  the  effect  of  station 
residual,  which  cat  introduce  large  uiases  in  the  estimation 
of  dT/dA.  Nor  are  the  effects  of  lateral  velocity  variations 
included  in  the  figures.  For  this  reason  these  theoretical 
results  should  be  regarded  as  lower  bounds  or  the  location 
uncertainty  to  be  attained  with  accurate  station  corrections. 

In  Figure  11  the  distribution  of  distance  and  azimuth 
errors  is  plotted  for  LASA  and  NORSAR.  All  earthquakes  in 
the  Summary  Bulletin  that  could  be  matched  with  ERL  epi¬ 
center  de terminations  for  the  specific  time  periods  were 
represented  in  the  distribution.  Events  with  large  distance 
and  azimuth  errors  were  generally  associated  with  events 
in  the  shadow  zone.  Biases  in  distance  and  azimuth  deter¬ 
minations  were  apparently  negligible  implying  that  station 
residual  effects  have  been  properly  compensated  for  in  the 
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Event  Processor.  The  magnitude  of  these  location  errors 
also  indicates  that  the  arrays  are  fairly  close  to  their 
theoretical  resolution  limit. 


DISCUSSION  AND  CONCLUSIONS 

In  this  study,  the  capabilities  of  LASA  and  NORSAR 
were  evaluated  on  the  basis  of  their  present  signal  pro¬ 
cessors.  The  effectiveness  of  this  analysis  was  limited 
by  several  factors.  First,  since  the  study  was  restricted 
to  data  that  had  been  partially  reduced  by  the  Detection 
Processor,  the  Event  Processor  and  the  analyst,  the  tendency 
would  be  to  underestimate  the  ultimate  capabilities  of  the 
arrays.  For  instance,  faulty  station  corrections  incorpor¬ 
ated  into  the  Detection  Processor  for  the  May  through  August 
1971  period  deteriorated  its  performance.  Second,  there  was 
a  lack  of  a  firm  and  large  standard  data  base  from  which  to 
evaluate  the  full  capabilities  of  LASA  and  NORSAR.  The 
only  set  of  known  earthquakes  available  at  the  time  of  the 
study  was  the  Environmental  Research  Laboratory  (ERL;  list 
of  epicenter  determinations.  This  list  omitted  many  low 
magnitude  events  which  LASA  and  NORSAR  were  capable  of 
observing.  Finally,  major  modifications  of  the  signal 
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processing  system  to  improve  its  performance  and  to  remove 
recently  discovered  errors  made  uniform  data  over  a  long 
time  span  unavailable. 

Notwithstanding  these  factors,  it  has  been  made 
apparent  that  excellent  detection  capability  can  be  achieved 
by  these  arrays,  LASA  in  particular.  The  location  capabili¬ 
ties  of  the  arrays,  although  not  as  good  as  those  of  a 
network  of  seismometers  adequately  spaced  around  the  epi¬ 
center,  are  close  to  their  theoretical  capability. 

The  determination  of  hypocenter  depth  by  a  single 
array  has  not  been  successful  due  to  the  problem  of  recog¬ 
nizing  and  distinguishing  the  depth  phases  pP  and  sP 
consistently.  For  this  reason  depth  determinations  were 
not  included  in  this  study. 

The  main  conclusions  of  the  study  are  listed  below. 

1.  LASA  detects  over  90  percent  of  the  events  re¬ 
ported  by  ERL  Ip  the  distance  range  30  to  SO 
degrees.  NORSAR  detects  about  60  percent  of  the 
reported  ERL  epicenter  determinations  in  the 
same  distance  range. 

2.  LASA's  detection  capability  is  at  least  0.5  magnitude 
lower  than  ERL's  capability  on  the  basis  of  the 
frequency  magnitude  distributions  of  the  LASA 
Summary  Bulletin  and  ERL  epicenter  determinations. 
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(ERL  reports  only  the  events  that  it  can  locate 
within  a  specified  accuracy.  There  are  many  more 
events  that  they  detect  but  do  not  list  in  their 
bulletin.)  NORSAR's  detection  capability  is  just 
slightly  better  than  ERL's. 

3.  By  extrapolation  of  the  Maximum  Short  Term  Average 
(MSTA)  of  detection  group  distribution  of  signals 
to  low  MSTA's  it  was  estimated  that  as  many  as  120 
seismic  signals  a  day  were  detected  by  LASA  and  35 
by  NORSAR.  The  problem  of  distinguishing  signals 
from  false  alarms  limits  the  number  of  seismic 
events  that  can  be  reported  to  those  with  high 
signal-to-noise  ratios. 

4.  LASA  and  NORSAR  are  able  to  locate  events  to  within 
3  degrees.  Errors  in  azimuth  determinations  are 
smaller  than  errors  in  distance  determination,  which 
is  in  agreement  with  the  theoretical  model  presented. 
It  appears  that  both  arrays  are  locating  events  at 
nearly  their  resolution  limit,  imposed  by  the  aper¬ 
ture  of  the  arrays. 

5.  NORSAR  underestimates  the  magnitude  of  an  earthquake 
by  one  half  a  unit  relative  to  ERL.  LASA's  magni¬ 
tude  bias  is  negligible  and  magnitude  dependent. 
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FIGURE  CAPTIONS 

.  Upper:  block  diagram  of  the  LASA  signal  processor. 

Lower:  block  diagram  of  the  LASA  Detection  Pro¬ 
cessor  . 

.  P  and  PKP  phase  locations  of  LASA  beams  (high 
resolution  partition)  as  given  by  the  Seismic 
Array  Analysis  Center  and  plotted  on  an  equi¬ 
distant  projection  centered  at  LASA. 

.  Upper:  cumulative  distribution  of  the  number  of 
detection  groups  (dots)  and  inferred  number  of 
signal  groups  (dashed  curve)  versus  MSTA  for  LASA 
(loft)  and  NORSAR  (right) ,  determined  for  the 
winter  period  1972.  Lower:  estimated  false  alarm 
probability  given  MSTA  for  LASA (left)  and  NORSAR 
(right).  Magnitude  given  for  earthquake  at  a  distance  of  60 

Percentage  of  ERL  events  detected  in  10  degree 
intervals  of  distance  for  LASA  (left)  and  NORSAR 
(right) .  Dots  and  solid  lines  denote  summer  1971 
period;  X's  and  dashed  lines  denote  winter  1972 
period. 

Percentage  of  ERL  events  detected  in  half  magnitude 
units  for  LASA  (left)  and  NORSAR  (right)  .  Upper 


Fig.  5. 
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plots  include  only  events  less  than  80  degrees 
distance  from  the  respective  array.  Lower 
plots  contain  the  events  further  than  80  degrees. 
Solid  line  -  summer  1971,  dashed  line  -  winter 
1972. 

Fig.  6.  Frequency  magnitude  distribution  of  earthquakes 

less  than  95  degrees  from  their  respective  arrays 
determined  from  LASA  and  NORSAR  Summary  Bulletins 
and  the  ERL  epicenter  determinations.  LASA's 
distribution  was  determined  for  summer  1971  and 
NORSAR' s  distribution  for  winter  1972. 

Fig.  7.  Correlation  plots  of  LASA  magnitude  versus  ERL 
magnitude  (summer  1971)  and  NORSAR  magnitude 
versus  ERL  magnitude  (winter  1972) . 

Fig.  8.  dT/dA  and  d^T/dA^  versus  A  for  P  phase. 

I ig .  9.  Distance  and  azimuth  of  epicenter  triggering 
specific  LASA  high  resolution  beams  assuming 
P  phase. 

Fig.  10.  Standard  error  in  azimuth  versus  a  . 

o 

Figl  11.  Histogram,  of  distance  and  azimuth  location  errors 
determined  for  LASA  (left  -  summer  1971)  and 
NORSAR  (right  -  winter  1972). 
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